Merge branch 'PhasicFlow:main' into stlModelImport

readme.md files update

.gitignore

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

CMakeLists.txt

@@ -1,75 +1,64 @@
 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 20 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 Release CACHE STRING "build type")
-
+set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
-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)
 
-option(USE_STD_PARALLEL_ALG "Use TTB std parallel algorithms" ON)
+message(STATUS "Install prefix is:" ${CMAKE_INSTALL_PREFIX})
+
+include(cmake/globals.cmake)
+
+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)
+#for installing the required packages
+include(cmake/preReq.cmake)
 
 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_DEFAULT_HOST_PARALLEL_EXECUTION_SPACE Serial CACHE STRING "" FORCE)
+	set(Kokkos_DEFAULT_DEVICE_PARALLEL_EXECUTION_SPACE OpenMP 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()
 
-
+if(pFlow_Build_MPI)
-
+	find_package(MPI REQUIRED)
-include(cmake/globals.cmake)
+endif()
-message(STATUS "Valid file extensions are ${validFiles}")
 
 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 +66,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(0)<<"Initializing host/device execution spaces . . . \n";
-	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
+	REPORT(1)<<"Host execution space is "<< Green_Text(DefaultHostExecutionSpace::name())<<END_REPORT;
-	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
+	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,273 @@
+/*------------------------------- 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)<<"Reading shape dictionary ..."<<END_REPORT;
+	grains_ = makeUnique<grainShape>(
+		pFlow::shapeFile__,
+		&Control().caseSetup(),
+		Property() );	
+
+	REPORT(0)<<"\nReading grain particles . . ."<<END_REPORT;
+	particles_ = makeUnique<grainFluidParticles>(Control(), grains_());
+
+
+	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,165 @@
+/*------------------------------- 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<grainShape> 			grains_   = 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 grainShape& grains)
+    : grainParticles(control, grains),
+      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() 
+{
+	const auto ti = this->TimeInfo();
+	accelerationTimer().start();
+		pFlow::grainFluidParticlesKernels::acceleration(
+			control().g(),
+			mass().deviceViewAll(),
+			contactForce().deviceViewAll(),
+			fluidForce_.deviceViewAll(),
+			I().deviceViewAll(),
+			contactTorque().deviceViewAll(),
+			fluidTorque_.deviceViewAll(),
+			pStruct().activePointsMaskDevice(),
+			acceleration().deviceViewAll(),
+			rAcceleration().deviceViewAll()
+			);
+	accelerationTimer().end();
+	
+	intCorrectTimer().start();
+	
+		dynPointStruct().correct(ti.dt());
+	
+		rVelIntegration().correct(ti.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 grainShape& grains);
+
+	~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,37 @@ 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)<<"Reading shapes dictionary..."<<END_REPORT;
+	spheres_ = makeUnique<sphereShape>(
+		pFlow::shapeFile__, 
+		&Control().caseSetup(), 
+		Property());
 
-	REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
+	REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-	particles_ = makeUnique<sphereFluidParticles>(Control(), Property());
+	particles_ = makeUnique<sphereFluidParticles>(Control(), spheres_());
 
 
 	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 +106,7 @@ pFlow::sphereDEMSystem::sphereDEMSystem(
 
 }
 
+
 pFlow::sphereDEMSystem::~sphereDEMSystem()
 {
 	
@@ -119,8 +128,8 @@ bool pFlow::sphereDEMSystem::updateParticleDistribution(
 	}
 
 	if(!particleDistribution_->locateParticles(
-			parPosition_, 
+			positionHost_, 
-			particles_->pStruct().activePointsMaskH()))
+			particles_->pStruct().activePointsMaskHost()))
 	{
 		fatalErrorInFunction<<
 		"error in locating particles among sub-domains"<<endl;
@@ -137,7 +146,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 +165,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 +232,14 @@ bool pFlow::sphereDEMSystem::sendFluidTorqueToDEM()
 
 bool pFlow::sphereDEMSystem::beforeIteration()
 {
-	parVelocity_ = particles_->dynPointStruct().velocityHostAll();
+	velocityHost_ = std::as_const(particles_()).velocity().hostView();
-	parPosition_ = particles_->dynPointStruct().pointPositionHostAll();
+	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
-	parDiameter_ = particles_->diameter().hostVectorAll();
+	diameterHost_ = particles_->diameter().hostView();
+
+	if(requireRVel_)
+		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();
+	
+
 	return true;
 }
 

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -42,24 +42,31 @@ protected:
 	
 	// protected members 
 
-	uniquePtr<property> 		property_ = nullptr;
+	uniquePtr<property> 			property_ = nullptr;
 
-	uniquePtr<geometry> 		geometry_ = nullptr;
+	uniquePtr<geometry> 			geometry_ = nullptr;
 
-	uniquePtr<sphereFluidParticles> 	particles_ = nullptr;
+	uniquePtr<sphereShape> 			spheres_  = nullptr;
 
-	uniquePtr<sphereInsertion> 	insertion_ = nullptr;
+	uniquePtr<sphereFluidParticles> particles_ = nullptr;
 
-	uniquePtr<interaction> 		interaction_ = nullptr;
+	uniquePtr<sphereInsertion> 		insertion_ = nullptr;
 
-	uniquePtr<domainDistribute> particleDistribution_=nullptr;
+	uniquePtr<interaction> 			interaction_ = 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 +99,8 @@ public:
 		word  demSystemName,
 		const std::vector<box>& domains,
 		int argc, 
-		char* argv[]);
+		char* argv[],
+		bool requireRVel=false);
 
 	virtual ~sphereDEMSystem();
 
@@ -110,15 +118,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,85 +21,85 @@ 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,
+    systemControl &control,
-	const property& prop
+    const sphereShape& shpShape)
-)
+    : sphereParticles(control, shpShape),
-:
+      fluidForce_(
-	sphereParticles(control, prop),
+          objectFile(
-	fluidForce_(
+              "fluidForce",
-		this->time().emplaceObject<realx3PointField_HD>(
+              "",
-			objectFile(
+              objectFile::READ_IF_PRESENT,
-				"fluidForce",
+              objectFile::WRITE_ALWAYS),
-				"",
+          dynPointStruct(),
-				objectFile::READ_IF_PRESENT,
+          zero3),
-				objectFile::WRITE_ALWAYS
+      fluidTorque_(
-				),
+          objectFile(
-			pStruct(),
+              "fluidTorque",
-			zero3
+              "",
-			)
+              objectFile::READ_IF_PRESENT,
-		),
+              objectFile::WRITE_NEVER),
-	fluidTorque_(
+          dynPointStruct(),
-		this->time().emplaceObject<realx3PointField_HD>(
+          zero3)
-			objectFile(
+{
-				"fluidTorque",
+	checkHostMemory();
-				"",
+} 
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			zero3
-			)
-		)
-{} 
 
 bool pFlow::sphereFluidParticles::beforeIteration()
 {
 	sphereParticles::beforeIteration();
+	checkHostMemory();
+
 	return true;
 }
 
 
 bool pFlow::sphereFluidParticles::iterate() 
 {
-
+	const auto ti = this->TimeInfo();
-	accelerationTimer_.start();
+	accelerationTimer().start();
 		pFlow::sphereFluidParticlesKernels::acceleration(
 			control().g(),
-			mass().deviceVectorAll(),
+			mass().deviceViewAll(),
-			contactForce().deviceVectorAll(),
+			contactForce().deviceViewAll(),
-			fluidForce().deviceVectorAll(),
+			fluidForce_.deviceViewAll(),
-			I().deviceVectorAll(),
+			I().deviceViewAll(),
-			contactTorque().deviceVectorAll(),
+			contactTorque().deviceViewAll(),
-			fluidTorque().deviceVectorAll(),
+			fluidTorque_.deviceViewAll(),
-			pStruct().activePointsMaskD(),
+			pStruct().activePointsMaskDevice(),
-			accelertion().deviceVectorAll(),
+			acceleration().deviceViewAll(),
-			rAcceleration().deviceVectorAll()
+			rAcceleration().deviceViewAll()
 			);
-	accelerationTimer_.end();
+	accelerationTimer().end();
 	
-	intCorrectTimer_.start();
+	intCorrectTimer().start();
 	
-		dynPointStruct_.correct(this->dt(), accelertion_);
+		dynPointStruct().correct(ti.dt());
 	
-		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
+		rVelIntegration().correct(ti.dt(), rVelocity(), rAcceleration());
 	
-	intCorrectTimer_.end();
+	intCorrectTimer().end();
 	
 	return true;
 }
 
 void pFlow::sphereFluidParticles::fluidForceHostUpdatedSync()
 {
-	fluidForce_.modifyOnHost();
+	copy(fluidForce_.deviceView(), fluidForceHost_);
-	fluidForce_.syncViews();
 	return;
 }
 
 void pFlow::sphereFluidParticles::fluidTorqueHostUpdatedSync()
 {
-	fluidTorque_.modifyOnHost();
+	copy(fluidTorque_.deviceView(), fluidTorqueHost_);
-	fluidTorque_.syncViews();
 	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,12 +61,12 @@ protected:
 	void zeroFluidTorque_H()
 	{
 		fluidTorque_.fillHost(zero3);
-	}
+	}*/
 
 public:
 
 	/// construct from systemControl and property 
-	sphereFluidParticles(systemControl &control, const property& prop);
+	sphereFluidParticles(systemControl &control, const sphereShape& shpShape);
 	
 	/// before iteration step 
 	bool beforeIteration() override;
@@ -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

@@ -1,40 +1,76 @@
-<div align ="center">
+<div align="center">
-<img src="doc/phasicFlow_logo_github.png" style="width: 400px;">
+  <img src="doc/phasicFlow_logo_github.png" style="width: 400px;" alt="PhasicFlow Logo">
 </div>
 
+## **PhasicFlow: High-Performance Discrete Element Method Simulations**
 
-**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. 
+PhasicFlow is a robust, open-source C++ framework designed for the efficient simulation of granular materials using the Discrete Element Method (DEM). Leveraging parallel computing paradigms, PhasicFlow is capable of executing simulations on shared-memory multi-core architectures, including CPUs and NVIDIA GPUs (CUDA-enabled). The core parallelization strategy focuses on loop-level parallelism, enabling significant performance gains on modern hardware. Users can seamlessly transition between serial execution on standard PCs, parallel execution on multi-core CPUs (OpenMP), and accelerated simulations on GPUs. Currently, PhasicFlow supports simulations involving up to 80 million particles on a single desktop workstation. Detailed performance benchmarks are available on the [PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Performance-of-phasicFlow).
 
-## How to build?
+**Scalable Parallelism: MPI Integration**
-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).
 
-## Online code documentation
+Ongoing development includes the integration of MPI-based parallelization with dynamic load balancing. This enhancement will extend PhasicFlow's capabilities to distributed memory environments, such as multi-GPU workstations and high-performance computing clusters. Upon completion, PhasicFlow will offer six distinct execution modes:
-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/)
 
-## How to use PhasicFlow?
+1.  **Serial Execution:** Single-core CPU.
-You can navigate into [tutorials folder](./tutorials) in the phasicFlow folder to see some simulation case setups. If you need more detailed discription, visit our [wiki page tutorials](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials).  
+2.  **Shared-Memory Parallelism:** Multi-core CPU (OpenMP).
+3.  **GPU Acceleration:** NVIDIA GPU (CUDA).
+4.  **Distributed-Memory Parallelism:** MPI.
+5.  **Hybrid Parallelism:** MPI + OpenMP.
+6.  **Multi-GPU Parallelism:** MPI + CUDA.
 
-## [PhasicFlowPlus](https://github.com/PhasicFlow/PhasicFlowPlus)
+## **Build and Installation**
-PhasicFlowPlus is and extension to PhasicFlow for simulating particle-fluid systems using resolved and unresolved CFD-DEM. [See the repository of this package.](https://github.com/PhasicFlow/PhasicFlowPlus)
 
+PhasicFlow can be compiled for both CPU and GPU execution.
 
-## Supporting packages
+* **Current Development (v-1.0):** Comprehensive build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-build-PhasicFlow%E2%80%90v%E2%80%901.0).
-* [Kokkos](https://github.com/kokkos/kokkos) from National Technology & Engineering Solutions of Sandia, LLC (NTESS)
+* **Latest Release (v-0.1):** Detailed build instructions are available [here](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-Build-PhasicFlow).
-* [CLI11 1.8](https://github.com/CLIUtils/CLI11) from University of Cincinnati.
+
+## **Comprehensive Documentation**
+
+In-depth documentation, including code structure, features, and usage guidelines, is accessible via the [online documentation portal](https://phasicflow.github.io/phasicFlow/).
+
+### **Tutorials and Examples**
+
+Practical examples and simulation setups are provided in the [tutorials directory](./tutorials). For detailed explanations and step-by-step guides, please refer to the [tutorial section on the PhasicFlow Wiki](https://github.com/PhasicFlow/phasicFlow/wiki/Tutorials).
+
+## Contributing to PhasicFlow
+We welcome contributions to PhasicFlow! Whether you're a developer or a new user, there are many ways to get involved. Here's how you can help:
+1. Bug Reports
+2. Suggestions for better user experience
+3. Feature request and algorithm improvements
+4. Tutorials, Simulation Case Setups and documentation
+5. Direct Code Contributions
+
+For more details on how you can contribute to PhasicFlow see [this page](https://github.com/PhasicFlow/phasicFlow/wiki/How-to-contribute-to-PhasicFlow). 
+
+## **PhasicFlowPlus: Coupled CFD-DEM Simulations**
+
+PhasicFlowPlus is an extension of PhasicFlow that facilitates the simulation of particle-fluid systems using resolved and unresolved CFD-DEM methods. The repository for PhasicFlowPlus can be found [here](https://github.com/PhasicFlow/PhasicFlowPlus).
+
+## How to cite PhasicFlow?
 
-## 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,
+@article
-title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
+{
-journal = {Computer Physics Communications},
+  NOROUZI2023108821,
-volume = {291},
+  title = {PhasicFlow: A parallel, multi-architecture open-source code for DEM simulations},
-pages = {108821},
+  journal = {Computer Physics Communications},
-year = {2023},
+  volume = {291},
-issn = {0010-4655},
+  pages = {108821},
-doi = {https://doi.org/10.1016/j.cpc.2023.108821},
+  year = {2023},
-url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
+  issn = {0010-4655},
-author = {H.R. Norouzi},
+  doi = {https://doi.org/10.1016/j.cpc.2023.108821},
-keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
+  url = {https://www.sciencedirect.com/science/article/pii/S0010465523001662},
+  author = {H.R. Norouzi},
+  keywords = {Discrete element method, Parallel computing, CUDA, GPU, OpenMP, Granular flow}
 }
 ```
+
+
+## **Dependencies**
+
+PhasicFlow relies on the following external libraries:
+
+* **Kokkos:** A community-led performance portability ecosystem within the Linux Foundation's High-Performance Software Foundation (HPSF). ([https://github.com/kokkos/kokkos](https://github.com/kokkos/kokkos))
+* **CLI11 1.8:** A command-line interface parser developed by the University of Cincinnati. ([https://github.com/CLIUtils/CLI11](https://github.com/CLIUtils/CLI11))

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.004);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.004;           // minimum space between centers of particles
+        
+        numPoints  1000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_1mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.006);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.006;           // minimum space between centers of particles
+        
+        numPoints  250000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_250KParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_2mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_2mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.003);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_2mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_2mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.2);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.2);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.2);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.2);         // second point
+        
+        p3        ( 0.2  0.2 1.2);         // third point
+        
+        p4        (-0.2  0.2 1.2);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  2000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.19);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_2mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.003);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_4mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_4mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.6);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.6);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.6);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.6);         // second point
+        
+        p3        ( 0.2  0.2 1.6);         // third point
+        
+        p4        (-0.2  0.2 1.6);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  4000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.59);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_4mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.005);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_500kParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_500kParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  0.8);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 0.8);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 0.8);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.8);         // second point
+        
+        p3        ( 0.2  0.2 0.8);         // third point
+        
+        p4        (-0.2  0.2 0.8);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.005;           // minimum space between centers of particles
+        
+        numPoints  500000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  0.79);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_500KParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/interaction

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     interaction;
+objectType     dictionary;
+fileFormat     ASCII;
+
+/*---------------------------------------------------------------------------*/
+
+materials   (glassMat wallMat); // a list of materials names
+
+densities   (2500.0 2500);      // density of materials [kg/m3]
+
+contactListType   sortedContactList;
+
+contactSearch
+{
+    method         NBS;
+
+    updateInterval   20;
+
+    sizeRatio       1.1;
+
+    cellExtent     0.55;
+
+    adjustableBox   Yes;
+}
+
+model
+{
+    contactForceModel    nonLinearLimited;
+
+    rollingFrictionModel normal;
+    
+    /*
+       Property (glassMat-glassMat      glassMat-wallMat
+                                        wallMat-wallMat);
+    */
+
+    Yeff (1.0e6  1.0e6  
+                 1.0e6);                                    // Young modulus [Pa]
+   
+    Geff (0.8e6  0.8e6  
+                 0.8e6);                                    // Shear modulus [Pa]
+
+    nu   (0.25   0.25  
+                 0.25);                                     // Poisson's ratio [-]
+
+    en   (0.97   0.85   
+                 1.00);                                     // coefficient of normal restitution
+
+    mu   (0.65   0.65   
+                 0.65);                                     // dynamic friction 
+
+    mur  (0.1    0.1    
+                 0.1);                                      // rolling friction 
+}
+

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/shapes

@@ -0,0 +1,15 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+
+names          (glassBead);    // names of shapes 
+
+diameters      (0.002);        // diameter of shapes 
+
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_8mParticles/cleanThisCase

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+ls | grep -P "^(([0-9]+\.?[0-9]*)|(\.[0-9]+))$" | xargs -d"\n" rm -rf
+rm -rf VTK 
+
+#------------------------------------------------------------------------------

benchmarks/rotatingDrum/rotatingDrum_8mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/domainDict

@@ -0,0 +1,50 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName     domainDict;
+objectType     dictionary;
+fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
+globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
+{
+    min (-0.2 -0.2  0.0);
+    max ( 0.2  0.2  1.6);
+}
+
+boundaries
+{
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
+    left
+    {
+        type    exit;                            // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;                           // other options: periodic, reflective 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/geometryDict

@@ -0,0 +1,86 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     geometryDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+motionModel rotatingAxis;         // motion model: rotating object around an axis
+
+rotatingAxisInfo                  // information for rotatingAxisMotion motion model 
+{
+    rotAxis 
+    {
+        p1    (0.0 0.0 0.0);      // first point for the axis of rotation    
+        
+        p2    (0.0 0.0 1.0);      // second point for the axis of rotation    
+        
+        omega 1.256;              // rotation speed (rad/s) => 12 rpm
+    }
+}
+
+surfaces
+{
+    cylinder
+    {
+        type        cylinderWall;          // type of the wall
+        
+        p1          (0.0 0.0 0.0);         // begin point of cylinder axis
+           
+        p2          (0.0 0.0 1.6);         // end point of cylinder axis
+        
+        radius1     0.2;                   // radius at p1
+        
+        radius2     0.2;                   // radius at p2
+        
+        resolution  24;                    // number of divisions
+        
+        material    wallMat;               // material name of this wall
+        
+        motion      rotAxis;               // motion component name 
+    }
+
+    /*
+        This is a plane wall at the rear end of cylinder
+    */
+    
+    wall1
+    {
+        type      planeWall;               // type of the wall  
+        
+        p1        (-0.2 -0.2 0.0);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 0.0);         // second point
+        
+        p3        ( 0.2  0.2 0.0);         // third point
+        
+        p4        (-0.2  0.2 0.0);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+
+    /*
+        This is a plane wall at the front end of cylinder
+    */
+    
+    wall2
+    {
+        type      planeWall;               // type of the wall
+        
+        p1        (-0.2 -0.2 1.6);         // first point of the wall
+        
+        p2        ( 0.2 -0.2 1.6);         // second point
+        
+        p3        ( 0.2  0.2 1.6);         // third point
+        
+        p4        (-0.2  0.2 1.6);         // fourth point 
+        
+        material  wallMat;                 // material name of the wall  
+        
+        motion    rotAxis;                 // motion component name 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/particlesDict

@@ -0,0 +1,47 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+
+objectName     particlesDict;
+objectType     dictionary;
+fileFormat     ASCII;
+
+setFields
+{
+    defaultValue 
+    {
+        velocity         realx3     (0 0 0);     // linear velocity (m/s)
+        acceleration     realx3     (0 0 0);     // linear acceleration (m/s2)
+        rotVelocity      realx3     (0 0 0);     // rotational velocity (rad/s)
+        shapeName        word       glassBead;   // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method ordered;
+
+    orderedInfo
+    {
+        diameter   0.003;           // minimum space between centers of particles
+        
+        numPoints  6000000;          // number of particles in the simulation 
+        
+        axisOrder  (z x y);         // axis order for filling the space with particles
+    }
+
+    regionType cylinder;            // other options: box and sphere    
+
+    cylinderInfo                    // cylinder for positioning particles 
+    {
+        p1     (0.0  0.0  0.01);    // lower corner point of the box 
+        
+        p2     (0.0  0.0  1.59);    // upper corner point of the box 
+        
+        radius 0.195;               // radius of cylinder 
+    }
+}

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/settingsDict

@@ -0,0 +1,34 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName    settingsDict;
+objectType    dictionary;
+fileFormat    ASCII;
+/*---------------------------------------------------------------------------*/
+run             rotatingDrum_4mParticles;
+
+dt              0.00001;      // time step for integration (s)
+
+startTime       0;            // start time for simulation 
+
+endTime         4;           // end time for simulation 
+
+saveInterval    0.2;          // time interval for saving the simulation
+
+timePrecision   5;            // maximum number of digits for time folder 
+
+g               (0 -9.8 0);   // gravity vector (m/s2) 
+
+includeObjects (diameter);                     // save necessary (i.e., required) data on disk
+
+// exclude unnecessary data from saving on disk
+excludeObjects (rVelocity.dy1 pStructPosition.dy1 pStructVelocity.dy1); 
+
+integrationMethod         AdamsBashforth2;     // integration method 
+
+writeFormat               binary;              // data writting format (ascii or binary)
+
+timersReport              Yes;      
+
+timersReportInterval      0.01;    

benchmarks/rotatingDrum_4MParticles/caseSetup/interaction

@@ -1,59 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName interaction;
-objectType dicrionary;
-
-materials   (glassMat wallMat); // a list of materials names
-densities   (2500.0 2500);      // density of materials [kg/m3]
-
-contactListType   sortedContactList;
-
-model
-{
-   contactForceModel nonLinearLimited;
-   rollingFrictionModel normal;
-
-   Yeff (1.0e6 1.0e6  // Young modulus [Pa]
-               1.0e6);
-   
-   Geff (0.8e6 0.8e6  // Shear modulus [Pa]
-               0.8e6);
-
-   nu    (0.25 0.25  // Poisson's ratio [-]
-               0.25);
-
-   en (0.97  0.85   // coefficient of normal restitution
-             1.00);
-
-   et (1.0   1.0    // coefficient of tangential restitution 
-             1.0);
-
-   mu (0.65  0.65   // dynamic friction 
-             0.65);
-
-   mur (0.1  0.1    // rolling friction 
-             0.1);
-        
-}
-
-contactSearch
-{
-   method         NBS;           
-   wallMapping    cellMapping;   
-
-   NBSInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list 
-      sizeRatio      1.05;       // bounding box size to particle diameter (max)
-   }
-
-   cellMappingInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list  
-      cellExtent     0.6;        // bounding box for particle-wall search (> 0.5)
-   }
-
-}

benchmarks/rotatingDrum_4MParticles/caseSetup/particleInsertion

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

benchmarks/rotatingDrum_4MParticles/caseSetup/sphereShape

@@ -1,11 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName 	sphereDict;
-objectType 	sphereShape;
-
-names 		(glassBead); 	// names of shapes 
-diameters 	(0.003);		// diameter of shapes 
-materials	(glassMat);		// material names for shapes 

benchmarks/rotatingDrum_4MParticles/settings/geometryDict

@@ -1,63 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	geometryDict;
-objectType 	dictionary;
-
-motionModel rotatingAxisMotion; 
-
-surfaces
-{
-	
-	cylinder
-	{
-		type cylinderWall;  
-		p1 (0.0 0.0 0.0);   
-		p2 (0.0 0.0 1.6);  
-		radius1 0.2;		
-		radius2 0.2;		
-		resolution 24;      
-		material wallMat;   
-		motion rotAxis;
-	}
-
-	
-	wall1
-	{
-		type planeWall;			
-		p1 (-0.2 -0.2 0.0);	
-		p2 ( 0.2 -0.2 0.0);   
-		p3 ( 0.2  0.2 0.0);   
-		p4 (-0.2  0.2 0.0);   
-		material wallMat;     
-		motion rotAxis;		
-	}
-
-	/*
-	This is a plane wall at the front end of cylinder
-	*/
-	wall2
-	{
-		type planeWall;
-		p1 (-0.2 -0.2 1.6);
-		p2 ( 0.2 -0.2 1.6);
-		p3 ( 0.2  0.2 1.6);
-		p4 (-0.2  0.2 1.6);
-		material wallMat;
-		motion rotAxis;
-	}
-		
-}
-
-// information for rotatingAxisMotion motion model 
-rotatingAxisMotionInfo
-{
-	rotAxis 
-	{
-		p1 (0.0 0.0 0.0);	
-		p2 (0.0 0.0 1.0);	
-		omega 1.256; 		// rotation speed (rad/s) => 12 rpm
-	}
-}

benchmarks/rotatingDrum_4MParticles/settings/particlesDict

@@ -1,44 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	particlesDict;
-objectType 	dictionary;
-
-setFields
-{
-	
-	defaultValue 
-	{
-		velocity 		realx3 	(0 0 0); // linear velocity (m/s)
-		acceleration 	realx3 	(0 0 0); // linear acceleration (m/s2)
-		rotVelocity 	realx3 	(0 0 0); // rotational velocity (rad/s)
-		shapeName 		word	glassBead; // name of the particle shape 
-	}
-
-	selectors
-	{}
-}
-
-positionParticles
-{
-	method positionOrdered;     		
-
-	maxNumberOfParticles 	4000001; 
-	mortonSorting 			Yes;    
-
-	cylinder  // box for positioning particles 
-	{
-		p1 ( 0.0  0.0 0.01);    // lower corner point of the box 
-		p2 ( 0.0  0.0 1.59);	// upper corner point of the box 
-		radius 0.195;
-	}
-	
-	positionOrderedInfo
-	{
-		diameter 0.003; 	// minimum space between centers of particles
-		numPoints 4000000; 	// number of particles in the simulation 
-		axisOrder (z x y);  // axis order for filling the space with particles
-	}
-}

benchmarks/rotatingDrum_4MParticles/settings/settingsDict

@@ -1,32 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-objectName 	settingsDict;
-objectType 	dictionary;;
-
-run 			rotatingDrum_1;
-
-dt 				0.00001; 	// time step for integration (s)
-
-startTime 		0; 			// start time for simulation 
-
-endTime 		10; 		// end time for simulation 
-
-saveInterval 	0.2; 		// time interval for saving the simulation
-
-timePrecision   5;			// maximum number of digits for time folder 
-
-g 				(0 -9.8 0); // gravity vector (m/s2) 
-
-domain 
-{
-	min (-0.2 -0.2 -0.0);
-	max ( 0.2  0.2  1.6);
-}
-
-integrationMethod 		AdamsBashforth3; 	// integration method 
-
-timersReport 			Yes;  	
-
-timersReportInterval   	0.01;	

cmake/autoComplete

@@ -0,0 +1,112 @@
+PF_cFlags="--description --help --version"
+AllTimeFolders=
+__getAllTime(){
+    # Initialize empty array for time folders
+    local time_folders=()
+    
+    # Loop through all directories in current folder
+    for dir in */; do
+        # Remove trailing slash
+        dir=${dir%/}
+        
+        # Check if directory name is a valid floating point number
+        # This pattern matches integers and floating point numbers
+        if [[ $dir =~ ^[0-9]+(\.[0-9]+)?$ ]]; then
+            time_folders+=("$dir")
+        fi
+    done
+    
+    # Set completion reply to the time folders
+    COMPREPLY=("${time_folders[@]}")
+    AllTimeFolders="${time_folders[@]}"
+}
+
+__getFields(){
+    __getAllTime
+    local -A unique_files=()
+    # Files to exclude from suggestions
+    local exclude_files=("shapeHash" "pStructure" "particleInsertion" "p" "alpha" "U" "Sp" "Su" "phi")
+    declare -A exclude_dict
+    
+    # Build exclude dictionary for faster lookups
+    for file in "${exclude_files[@]}"; do
+        exclude_dict["$file"]=1
+    done
+
+    for dir in $AllTimeFolders; do
+        # Skip if not a directory
+        [ ! -d "$dir" ] && continue
+        
+        # Find all files in this time directory
+        while IFS= read -r filename; do
+            # Skip empty lines and excluded files
+            [ -z "$filename" ] || [ "${exclude_dict[$filename]+exists}" ] && continue
+            
+            # Add to unique files
+            unique_files["$filename"]=1
+        done < <(find "$dir" -maxdepth 1 -type f -printf '%f\n')
+    done
+    
+    # Set completion reply to the unique filenames
+    COMPREPLY=(${!unique_files[@]})
+    
+    # Clear global variable
+    AllTimeFolders=
+}
+
+_pFlowToVTK(){
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    
+    # Check if we're completing a field
+    local is_field=0
+    for ((i=1; i<COMP_CWORD; i++)); do
+        if [[ "${COMP_WORDS[i]}" == "--fields" ]]; then
+            is_field=1
+            break
+        fi
+    done
+    
+    if [ "$prev" == "--time" ]; then
+        __getAllTime
+    elif [ "$prev" == "--fields" ] || [ $is_field -eq 1 ]; then 
+        # We're completing field names
+        __getFields
+        # Filter the results based on the current word prefix
+        if [ -n "$cur" ]; then
+            local filtered=()
+            for item in "${COMPREPLY[@]}"; do
+                if [[ "$item" == "$cur"* ]]; then
+                    filtered+=("$item")
+                fi
+            done
+            COMPREPLY=("${filtered[@]}")
+        fi
+    else
+        COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$cur") )
+    fi
+}
+complete -F _pFlowToVTK pFlowToVTK
+
+_postprocessPhasicFlow(){
+    if [ "$3" == "--time" ]; then
+      __getAllTime
+    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_VERSION=v-1.0
-
-export pFlow_PROJECT=phasicFlow
 
+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} source files are: ${${source_files}}")
-message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n")
+#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} source files are: ${source_files}")
-message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n")
+#message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n \n")
 
 install(TARGETS ${target_name} DESTINATION lib)
 install(FILES ${includeFiles}  DESTINATION include/${target_name})

cmake/preReq.cmake

@@ -0,0 +1,54 @@
+if(pFlow_STD_Parallel_Alg)
+    # Check if libtbb-dev is installed
+    execute_process(
+        COMMAND dpkg -s libtbb-dev
+        RESULT_VARIABLE TBB_IS_INSTALLED
+        OUTPUT_QUIET
+        ERROR_QUIET)
+
+    if(NOT TBB_IS_INSTALLED EQUAL 0)
+        message(STATUS "libtbb-dev not found. Installing libtbb-dev...")
+        execute_process(
+            COMMAND sudo apt-get update
+            COMMAND sudo apt-get install -y libtbb-dev
+            RESULT_VARIABLE TBB_INSTALL_RESULT)
+            
+        if(NOT TBB_INSTALL_RESULT EQUAL 0)
+            message(FATAL_ERROR "Failed to install libtbb-dev")
+        endif()
+    else()
+        message(STATUS "libtbb-dev is already installed.")
+    endif()
+endif()
+
+
+# Kokkos folder creation 
+set(Kokkos_Source_DIR $ENV{HOME}/Kokkos/kokkos)
+
+if(NOT EXISTS "${Kokkos_Source_DIR}/CMakeLists.txt")
+
+    # Check CMake version and set policy CMP0169 if CMake version is 3.30 or higher
+    if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
+        cmake_policy(SET CMP0169 OLD)
+    endif()
+
+    include(FetchContent)
+    FetchContent_Declare(
+        kokkos
+        GIT_REPOSITORY https://github.com/kokkos/kokkos.git
+        GIT_TAG 4.4.01
+    )
+    
+    FetchContent_GetProperties(kokkos)
+    if(NOT kokkos_POPULATED)
+        message(STATUS "Kokkos source directory not found. Downloading Kokkos version 4.4.1 ...")
+        FetchContent_Populate(kokkos)
+        set(Kokkos_Source_DIR ${kokkos_SOURCE_DIR})
+    endif()
+
+endif()
+
+message(STATUS "Kokkos source directory is ${Kokkos_Source_DIR}")
+add_subdirectory(${Kokkos_Source_DIR} ./kokkos)
+#Kokkos_cmake_settings()
+

doc/mdDocs/codingStyle.md

@@ -0,0 +1,231 @@
+# PhasicFlow Coding Style Guidelines
+This document outlines the coding style guidelines for the PhasicFlow codebase. 
+Adhering to these guidelines ensures consistency, readability, and maintainability of the project.
+## 1. FormattingIndentation: 
+* Use 4 spaces for every logical level and block.
+* Line Spacing: Leave two empty lines between sections (e.g., between functions in a .cpp file, between class members).
+
+## 2. Naming ConventionsGeneral Naming: 
+* All names should start with lowercase letters, except for special names (e.g., Ergun, Hertz).
+* Macro names start with Upper case or all the letters are in UPPER case.
+* Compound Names: For compound names, the first word starts with a lowercase letter, and subsequent words start with an uppercase letter (e.g., boundaryBase, motionModel).
+
+## 3. File Structure 
+* Header Files: Use the .hpp extension for header files.
+* Source Files: Use the .cpp extension for source files.
+* Header and Source File Headers: All header and source files must include a standardized header that describes the project's intention and licensing information.
+* File Naming: Header and source file names should correspond to the class they contain. Aim for one class per file.
+* Inline Functions: Place inline functions in a separate classNameI.hpp file to avoid cluttering the main header file.
+
+## 4. Class DesignClass Member Order:
+* Private members and methods
+* Private static members and methods
+* Public methods
+* Public static methods
+* Enumerations and Nested Classes: Declare enumerations and nested classes before all class members and methods.
+* Special Functions: Each class must explicitly define all special functions:Constructor, Copy constructor and assignment operator, Move constructor and assignment operator
+* Destructor: Each class must have an explicit destructor declaration:`~className() = default`; or `~className();`
+* Interface classes or classes with virtual methods must have a virtual destructor.
+* Virtual Method Overrides: When implementing a `virtual` method from a base class in a derived class, use the `override` keyword. The same applies to derived class destructors.
+## 5. NamespacesOfficial Namespace: 
+The official namespace for the codebase is pFlow. All entities should be defined within this namespace.
+
+### Example File Structure
+```
+src/
+├── componentName1/
+│   ├── componentName1.hpp
+│   ├── componentName1.cpp
+│   ├── componentName1I.hpp
+│   └── ...
+└── componentName2/
+    ├── componentName2.hpp
+    ├── componentName2.cpp
+    └── ...
+```
+### Example Class Structure
+```C++
+namespace pFlow 
+{
+
+class MyClass 
+{
+public:
+    enum class MyEnum 
+    {
+        Value1,
+        Value2
+    };
+
+    class NestedClass 
+    {
+        // ...
+    };
+
+private:
+    int privateMember_;
+    
+    void privateMethod();
+
+    
+    static int privateStaticMember_;
+    
+    static void privateStaticMethod();
+
+public:
+    
+    MyClass();
+
+    MyClass(const MyClass& other);
+    
+    MyClass(MyClass&& other);
+    
+    MyClass& operator=(const MyClass& other);
+    
+    MyClass& operator=(MyClass&& other);
+    
+    ~MyClass();
+
+    
+    void publicMethod();
+    
+    static void publicStaticMethod();
+
+};
+
+// assuming base class has virtual methods 
+class DerivedClass 
+: 
+    public BaseClass 
+{
+public:
+     
+     ...
+
+    ~DerivedClass() override;
+    
+     void virtualMethod() override;
+};
+
+} // namespace pFlow
+``` 
+
+## 6. Doxygen Documentation
+
+### 6.1. Ruls
+provide the documentations in the header files only. In rare cases where you are generating documentations for executables, the doxygen documentation can be supplied in the .cpp file.  
+
+* **General Doxygen Style:**
+
+    * Use `///` for short documentations for methods and members.
+
+    * Use `/** */` for classes and main methods which play a significant role in the class or code.
+
+    * Place Doxygen comments *before* the declaration of the entity being documented (e.g., class, function, variable).
+
+    * Use `@param` to document function parameters, `@return` for return values, `@brief` for a short description, and `@details` for a more in-depth explanation.
+
+    * Use Markdown syntax within Doxygen comments for formatting.
+
+* **File Headers:** Each file should contain a Doxygen comment at the top, including:
+
+    * `@file` : The name of the file.
+
+    * `@brief`: A brief description of the file's purpose.
+
+    * `@author`: The author(s) of the file.
+
+    * `@date` : The date of creation or last modification.
+
+* **Class Documentation:**
+
+    * Use `/** */` for class documentation.
+
+    * Provide a `@brief` description of the class.
+
+    * Use `@tparam` to document template parameters.
+
+    * Document the purpose of the class, its invariants, and how it should be used.
+
+* **Function/Method Documentation:**
+
+    * Use `///` for short documentations.
+
+    * Use `/** */` for main methods which play a significant role.
+
+    * Provide a `@brief` description of the function.
+
+    * Use `@param` to describe each parameter, including its purpose and whether it is an input, output, or input/output parameter.
+
+    * Use `@return` to describe the return value, including its meaning and possible values.
+
+    * Use `@pre` to document any preconditions that must be met before calling the function.
+
+    * Use `@post` to document any postconditions that will be true after the function returns.
+
+    * Use `@throws` to document any exceptions that the function may throw.
+
+    * Use `@details` for a more detailed explanation of the function's behavior, algorithms, or any other relevant information.
+
+* **Variable Documentation:**
+
+    * Use `///<` for single-line documentation of variables.
+
+### 6.2. Doxygen Documentation Examples
+
+* **Class example**
+
+```cpp
+/**
+ * @brief Represents a particle in the simulation.
+ * @details This class stores the position, velocity, and other physical
+ * properties of a particle.
+ */
+class Particle
+{
+private:
+
+    Point     position_;   ///< The current position of the particle.
+
+    Vector    velocity_;   ///< The current velocity of the particle.
+
+    double     mass_;      ///< The mass of the particle.
+
+public:
+
+    /// Constructs a particle with default values.
+    Particle();
+
+    /**
+     * @brief Updates the position of the particle based on its velocity
+     * and the given time step.
+     * @param deltaTime The time elapsed since the last update, in seconds.
+     */
+    void updatePosition(const timeInfo& ti );
+
+    /// Gets the current position of the particle.
+    Point getPosition() const;
+};
+```
+
+* **Function Example**
+
+```cpp
+/**
+ * @brief Calculates the distance between two points.
+ * @param p1 The first point.
+ * @param p2 The second point.
+ * @return The distance between the two points.
+ */
+double calculateDistance(const Point& p1, const Point& p2)
+{
+    // Implementation
+    return 0.0;
+}
+
+/// Returns the velocity of the particle.
+Vector getVelocity() const
+{
+    return velocity_;
+}
+``` 

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

@@ -17,35 +17,41 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
+REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
-template<typename T, typename... properties>
+pFlow::grainShape grains
-void insertSetElementH
 (
-	ViewType1D<T, properties...>& view,
+	pFlow::shapeFile__,
-	hostViewType1D<label>& selected,
+	&Control.caseSetup(),
-	T val
+	proprties
 );
 
-template<typename T, typename... properties>
+//
-void insertSetElementH
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-(
+pFlow::grainParticles grnParticles(Control, grains);
-	ViewType1D<T, properties...>& view,
-	hostViewType1D<label>& selected,
-	hostViewType1D<T>& vals
-);
 
-template<typename T, typename... properties>
+//
-void insertSetElementD
+REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
-(
+/*auto& sphInsertion = 
-	ViewType1D<T, properties...>& view,
+	Control.caseSetup().emplaceObject<sphereInsertion>(
-	deviceViewType1D<label>& selected,
+		objectFile(
-	T val
+			insertionFile__,
-);
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
+			),
+		sphParticles,
+		sphParticles.shapes()
+	);*/
 
-template<typename T, typename... properties>
+auto grnInsertion =  pFlow::grainInsertion(
-void insertSetElementD
+	grnParticles, 
-(
+	grnParticles.grains());
-	ViewType1D<T, properties...>& view,
+
-	deviceViewType1D<label>& selected,
+REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
-	deviceViewType1D<T>& vals
+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"
+//#include "readControlDict.hpp"
-
-using pFlow::commandLine;
 
 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

@@ -0,0 +1,33 @@
+/*------------------------------- 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)<<"Reading shape dictionary ..."<<END_REPORT;
+pFlow::sphereShape spheres
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
+
+//
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
+pFlow::sphereParticles sphParticles(Control, spheres);
+
+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

@@ -17,13 +17,21 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
+REPORT(0)<<"Reading shapes dictionary..."<<END_REPORT;
+pFlow::sphereShape spheres
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
+//
+REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
+
+
+pFlow::sphereParticles sphParticles(Control, spheres);
 
 //
-REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
+REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
-sphereParticles sphParticles(Control, proprties);
-
-//
-REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
 /*auto& sphInsertion = 
 	Control.caseSetup().emplaceObject<sphereInsertion>(
 		objectFile(
@@ -36,13 +44,12 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
 		sphParticles.shapes()
 	);*/
 
-auto sphInsertion =  sphereInsertion(
+auto sphInsertion =  pFlow::sphereInsertion(
-	Control.caseSetup().path()+insertionFile__, 
 	sphParticles, 
-	sphParticles.shapes());
+	sphParticles.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;
-auto interactionPtr = interaction::create(
+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 
 	{
-		
+		//Ping;
 		if(! sphInsertion.insertParticles( 
+			Control.time().currentIter(),
 			Control.time().currentTime(),
 			Control.time().dt()	) )
 		{
@@ -95,30 +85,38 @@ if(!cmds.parse(argc, argv)) return 0;
 			return 1;
 		}
 		
+		// set force to zero
 		surfGeometry.beforeIteration();
 
-		sphInteraction.beforeIteration();
+		// set force to zero, predict, particle deletion and etc. 
-
 		sphParticles.beforeIteration();
-		
+		//Ping;
+		sphInteraction.beforeIteration();
 		
 		sphInteraction.iterate();
 			
-		sphParticles.iterate();
-
 		surfGeometry.iterate();
 
-		sphParticles.afterIteration();
+		//Ping;
+		sphParticles.iterate();
 		
+		//Ping;
+		sphInteraction.afterIteration();
+
+		//Ping;
 		surfGeometry.afterIteration();
 
+		//Ping;
+		sphParticles.afterIteration();
+		//Ping;	
 
 	}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,11 @@ add_subdirectory(Property)
 
 add_subdirectory(Particles)
 
+add_subdirectory(Geometry)
+
 add_subdirectory(Interaction)
 
 add_subdirectory(MotionModel)
 
-add_subdirectory(Geometry)
+add_subdirectory(PostprocessData)
 

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()
 {
-
+	if(materialName_.size() != numSurfaces() )
-	int8Vector propId(0, surface().capacity(),RESERVE());
-	propId.clear();
-
-	uint32 pId;
-	ForAll(matI, materialName_)
 	{
+		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<<
+			fatalErrorInFunction<<
-		 	"material name for the geometry is invalid: "<< materialName_[matI]<<endl;
+			"Property/material name is invalid "<<materialName_[i]<<
-		 	return false;
+			". 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()),
+	propertyId_
-	triSurface_(
+	(
-		control.geometry().emplaceObject<multiTriSurface>(
+		objectFile
-			objectFile(
+		(
-				"triSurface",
+			"propertyId",
-				"",
+			"",
-				objectFile::READ_ALWAYS,
+			objectFile::READ_NEVER,
-				objectFile::WRITE_ALWAYS
+			objectFile::WRITE_NEVER
-				)
-			)
 		),
-	motionComponentName_(
+		*this,
-		control.geometry().emplaceObject<wordField>(
+		0u
-			objectFile(
+	),
-				"motionComponentName",
+	contactForceWall_
-				"",
+	(
-				objectFile::READ_ALWAYS,
+		objectFile
-				objectFile::WRITE_ALWAYS	
+		(
-				),
+			"contactForcWall",
-			"motionNamesList"
+			"",
-			)
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
 		),
-	materialName_(
+		*this,
-		control.geometry().emplaceObject<wordField>(
+		zero3
-			objectFile(
+	),
-				"materialName",
+	normalStressWall_
-				"",
+	(
-				objectFile::READ_ALWAYS,
+		objectFile
-				objectFile::WRITE_ALWAYS
+		(
-				),
+			"normalStressWall",
-			"materialNamesList"
+			"",
-			)
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_NEVER
 		),
-	propertyId_(
+		*this,
-		control.geometry().emplaceObject<int8TriSurfaceField_D>(
+		zero3
-			objectFile(
+	),
-				"propertyId",
+	shearStressWall_
-				"",
+	(
-				objectFile::READ_NEVER,
+		objectFile
-				objectFile::WRITE_NEVER),
+		(
-			surface(),
+			"shearStressWall",
-			0 ) ),
+			"",
-	contactForceWall_(
+			objectFile::READ_IF_PRESENT,
-		control.geometry().emplaceObject<realx3TriSurfaceField_D>(
+			objectFile::WRITE_NEVER
-			objectFile(
+		),
-				"contactForceWall",
+		*this,
-				"",
+		zero3
-				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) )
 {
 
-	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 
+		/// Stress on each triangle in the set of wall surfaces 
-		realx3TriSurfaceField_D&    stressWall_;
+		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()
+		void zeroForce();
-		{
-			contactForceWall_.fill(zero3);
-		}
-
 		
 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,
+				systemControl& control, 
-			 const dictionary& dict,
+				const property& prop,
-			 const multiTriSurface& triSurface,
+				multiTriSurface& surf,
-			 const wordVector& motionCompName,
+				const wordVector& motionCompName,
-			 const wordVector& propName),
+				const wordVector& materialName,
-			(control, prop, dict, triSurface, motionCompName, propName)
+				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();
+			return motionComponentName_;
-		}
-
-		/// 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_;
 		}
 		
 		/// 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
+		bool write( iOstream& os, const IOPattern& iop )const override;
-		{
+		
-			return owner().write();
-		}
 
 
 	//- 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,
+			systemControl& control, 
-				const property& prop,
+			const property& prop,
-				const dictionary& dict,
+			multiTriSurface& surf,
-				const multiTriSurface& triSurface,
+			const wordVector& motionCompName,
-				const wordVector& motionCompName,
+			const wordVector& materialName,
-				const wordVector& propName);
+			const dictionary& motionDic);
 
 };
 

src/Geometry/geometryMotion/geometryMotion.cpp

@@ -21,41 +21,115 @@ Licence:
 template<typename MotionModel>
 bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 {
-	motionIndex_.clear();
+	if(motionComponentName().size() != numSurfaces() )
-	triMotionIndex_.reserve( this->surface().capacity() );
-	triMotionIndex_.clear();
-
-	ForAll( surfI, motionComponentName_)
 	{
-		auto mName = motionComponentName_[surfI];
+		fatalErrorInFunction<<
-		auto mInd = motionModel_.nameToIndex(mName);
+		"size of motion component names in the triSurface is not"<<
-		motionIndex_.push_back(mInd);
+		" equal to size of number of sub-surfaces"<<endl;
-		// fill motionIndex for triangles of the surface
+		 return false;
-		int32 surfSize = this->surface().surfNumTriangles(surfI);
+	}
-		for(int32 i=0; i<surfSize; i++)
+
+	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());
+	surfMotionIndex_.assign(surfMotionIndex);
-	pointMotionIndex_.clear();
+	triMotionIndex_.assign(triMotionIndex);
-
+	pointMotionIndex_.assign(pointMotionIndex);
-	ForAll(surfI, motionIndex_)
-	{
-		auto nP = triSurface_.surfNumPoints(surfI);
-		for(int32 i=0; i<nP; i++)
-		{
-			pointMotionIndex_.push_back(motionIndex_[surfI]);
-		}
-	}
-	pointMotionIndex_.syncViews();
 
 	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 +139,75 @@ pFlow::geometryMotion<MotionModel>::geometryMotion
 )
 :
 	geometry(control, prop),
-	motionModel_(
+	motionModel_
-		this->owner().template emplaceObject<MotionModel>(
+	(
-			objectFile(
+		objectFile
-				motionModelFile__,
+		(
-				"",
+			motionModelFile__,
-				objectFile::READ_ALWAYS,
+			"",
-				objectFile::WRITE_ALWAYS
+			objectFile::READ_ALWAYS,
-				)
+			objectFile::WRITE_ALWAYS
-			)
 		),
+		owner()
+	),
 	moveGeomTimer_("move geometry", &this->timers())
 {
-	findMotionIndex();
+	if(!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() )
 	{
-		moveGeomTimer_.start();
+		fatalExit;
-		moveGeometry();
-		moveGeomTimer_.end();
 	}
-	return true;
 }
 
-template<typename MotionModel>
+template <typename MotionModelType>
-bool pFlow::geometryMotion<MotionModel>::afterIteration() 
+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();
+	if(!findMotionIndex())
-	return true;
+	{
+		fatalExit;
+	}
 } 
 
-template<typename MotionModel>
-bool pFlow::geometryMotion<MotionModel>::moveGeometry()
-{
 
-	real dt = this->dt();
+ template<typename MotionModel>
-	real t = this->currentTime();
+ bool pFlow::geometryMotion<MotionModel>::iterate()
-
+ {
-	auto pointMIndex= pointMotionIndex_.deviceVector();
+     if( motionModel_.isMoving() )
-	auto mModel = motionModel_.getModel(t);
+     {
-	realx3* points = triSurface_.pointsData_D();
+        moveGeomTimer_.start();
-	auto numPoints = triSurface_.numPoints();
+          moveGeometry();
-
+          this->calculateNormals();
-
+        moveGeomTimer_.end();
-	Kokkos::parallel_for(
+     }
-		"geometryMotion<MotionModel>::movePoints",
+     return true;
-		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;
-}

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 wordVector& materialName,
-			const MotionModel& motionModel);
+			const dictionary& motionDict);
-
-		/// 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);
 
 		/// 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"
+template class pFlow::geometryMotion<pFlow::vibratingMotion>;
-#endif
+
+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,8 +22,9 @@ Licence:
 #define __geometryMotions_hpp__
 
 #include "geometryMotion.hpp"
-#include "fixedWall.hpp"
+#include "stationaryWall.hpp"
 #include "rotatingAxisMotion.hpp"
+#include "conveyorBeltMotion.hpp"
 #include "multiRotatingAxisMotion.hpp"
 #include "vibratingMotion.hpp"
 
@@ -31,15 +32,15 @@ Licence:
 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>;
-
-typedef geometryMotion<fixedWall> fixedGeometry;
 
+using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
 
+using multiRotationAxisMotionGeometry = geometryMotion<multiRotatingAxisMotion>;
 
 }
 

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,87 +19,192 @@ 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
-
-pFlow::AdamsBashforth2::AdamsBashforth2
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	dy1_(
-		owner.emplaceObject<realx3PointField_D>(
-			objectFile(
-				groupNames(baseName,"dy1"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			zero3))
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth2::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth2::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
-{
+	real damping = 1.0)
-	if(this->pStruct().allActive())
-	{
-		return intAll(dt, y, dy, this->pStruct().activeRange());
-	}
-	else
-	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth2::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth2::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
 {
 
-	auto d_dy = dy.deviceVectorAll();
+	auto d_dy = dy.deviceView();
-	auto d_y  = y.deviceVectorAll();
+	auto d_y  = y.deviceView();
-	auto d_dy1= dy1_.deviceVectorAll();
+	auto d_dy1= dy1.deviceView();
+	auto activeRng = dy1.activeRange();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth2::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){
+		LAMBDA_HD(uint32 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_y[i] +=  damping * 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,
+	real damping = 1.0
+)
+{
+
+	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] +=  damping * 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,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField,
+	bool  keepHistory
+)
+:
+	integration(baseName, pStruct, method, initialValField, keepHistory),
+	realx3PointField_D
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy1"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
+		),
+		pStruct,
+		zero3,
+		zero3
+	),
+	initialValField_(initialValField),
+	boundaryList_(pStruct, method, *this)
+{
+	realx3PointField_D::addEvent(message::ITEMS_INSERT);
+}
+
+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 dt, 
+	realx3Field_D &y, 
+	realx3PointField_D &dy
+)
+{
+    return true;
+}
+
+bool pFlow::AdamsBashforth2::correct
+(
+	real dt,
+	realx3PointField_D& y,
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	auto& dy1l = dy1();
+	bool success = false;
+	if(dy1l.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), damping);
+	}
+
+	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,  dy1());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel,  dy1());
+	}
+
+	success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth2::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		
+		return insertValues(varList, initialValField_.deviceViewAll(), dy1());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+    
+}*/

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,56 @@ namespace pFlow
  */
 class AdamsBashforth2
 :
-	public integration
+	public integration,
+	public realx3PointField_D
 {
+private:
+
+	const realx3Field_D& initialValField_;
+
+	boundaryIntegrationList boundaryList_;
+
+	friend class processorAB2BoundaryIntegration;
+
 protected:
 	
-	/// dy at t-dt
+	const auto& dy1()const
-	realx3PointField_D& dy1_;
+	{
+		return static_cast<const realx3PointField_D&>(*this);
+	}
 
-	/// Range policy for integration kernel (alias)
+	auto& dy1()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return static_cast<realx3PointField_D&>(*this);
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
+
-			>;
+	auto& initialValField()
+	{
+		return initialValField_;
+	}
+
+	boundaryIntegrationList& boundaryList()
+	{
+		return boundaryList_;
+	}
 	
 public:
 
-	/// Type info
+	/// Class info
-	TypeInfo("AdamsBashforth2");
+	ClassInfo("AdamsBashforth2");
 
 	// - Constructors
 
 		/// Construct from components
 		AdamsBashforth2(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField,
-
+			bool  keepHistory);
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth2>(*this);
-		}
 
 		/// Destructor 
-		virtual ~AdamsBashforth2()=default;
+		~AdamsBashforth2()override = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -81,70 +96,44 @@ public:
 
 	// - Methods
 
+		void updateBoundariesSlaveToMasterIfRequested()override;
+
+		/// return integration method 
+		word method()const override
+		{
+			return "AdamsBashforth2";
+		}
+
 		bool predict(
 			real UNUSED(dt), 
-			realx3Vector_D& UNUSED(y), 
+			realx3PointField_D& UNUSED(y), 
-			realx3Vector_D& UNUSED(dy)) override;
+			realx3PointField_D& UNUSED(dy)) final;
+		
+		bool predict(
+			real dt, 
+			realx3Field_D& y, 
+			realx3PointField_D& dy) final;
 
 		bool correct(
 			real dt, 
-			realx3Vector_D& y, 
+			realx3PointField_D& y, 
-			realx3Vector_D& dy) override;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
-		bool setInitialVals(
+		bool correctPStruct(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
-		{
-			return false;
-		}
-		
-		/// Integrate on all points in the active range 
-		bool intAll(
 			real dt, 
-			realx3Vector_D& y, 
+			pointStructure& pStruct, 
-			realx3Vector_D& dy, 
+			realx3PointField_D& vel) override;
-			range activeRng);
 		
-		/// Integrate on active points in the active range 
+		/*bool hearChanges
-		template<typename activeFunctor>
+		(
-		bool intRange(
+			const timeInfo& ti,
-			real dt, 
+			const message& msg, 
-			realx3Vector_D& y, 
+			const anyList& varList
-			realx3Vector_D& dy, 
+		) override;*/
-			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