global damping is activated for velocity and rVelocity in both sphere and grain solvers

merge from main

CMakeLists.txt

@@ -6,7 +6,7 @@ project(phasicFlow VERSION 1.0 )
 set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
-set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
+set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE)
 set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
 mark_as_advanced(FORCE var BUILD_SHARED_LIBS)
 
@@ -23,7 +23,7 @@ else()
   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(${Kokkos_Source_DIR} ./kokkos)
 Kokkos_cmake_settings()
 
 
@@ -68,14 +68,13 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 #add a global include directory 
 include_directories(src/setHelpers src/demComponent "${PROJECT_BINARY_DIR}")
 
-
 add_subdirectory(src)
 
 add_subdirectory(solvers)
 
 add_subdirectory(utilities)
 
-#add_subdirectory(DEMSystems)
+add_subdirectory(DEMSystems)
 
 install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"
   DESTINATION include)

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",

cmake/autoComplete

@@ -0,0 +1,71 @@
+PF_cFlags="--description --help --version"
+AllTimeFolders=
+__getAllTime(){
+	files=( $(ls) )
+	deleteFiles=(settings caseSetup cleanThisCase VTK runThisCase stl postprocess postProcess) 
+	declare -A delk
+	for del in "${deleteFiles[@]}" ; do delk[$del]=1 ; done
+		# Tag items to remove, based on
+	for k in "${!files[@]}" ; do
+		[ "${delk[${files[$k]}]-}" ] && unset 'files[k]'
+	done
+		# Compaction
+	COMPREPLY=("${files[@]}")
+	AllTimeFolders="${files[@]}"
+}
+
+__getFields(){
+	__getAllTime
+	local -A unique_files=()
+
+	for dir in $AllTimeFolders; do
+	  # Check if the directory exists
+	  if [ ! -d "$dir" ]; then
+		continue  # Skip to the next directory
+	  fi
+		  
+	  files_in_dir=$(find "$dir" -maxdepth 1 -type f -printf '%f\n' | sort -u)
+	  
+	# Add filenames to the associative array (duplicates will be overwritten)
+	  while IFS= read -r filename; do
+		unique_files["$filename"]=1 # Just the key is important, value can be anything
+	  done <<< "$files_in_dir"
+	  
+	done
+	COMPREPLY=("${!unique_files[@]}")
+	AllTimeFolders=
+}
+
+_pFlowToVTK(){
+    if [ "$3" == "--time" ]  ; then
+        __getAllTime
+    elif [ "$3" == "--fields" ]; then 
+    	__getFields
+    else
+        COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") )
+    fi
+}
+complete -F _pFlowToVTK pFlowToVTK
+
+_postprocessPhasicFlow(){
+    if [ "$3" == "--time" ]; then
+      __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

@@ -3,7 +3,6 @@ export pFlow_PROJECT_VERSION=v-1.0
 
 export pFlow_PROJECT="phasicFlow-$pFlow_PROJECT_VERSION"
 
-
 projectDir="$HOME/PhasicFlow" 
 kokkosDir="$HOME/Kokkos/kokkos"
 
@@ -31,5 +30,5 @@ 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/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 \n")
+#message(STATUS "  ${target_name} include dirs are: ${includeDirs}\n \n")
 
 install(TARGETS ${target_name} DESTINATION lib)
 install(FILES ${includeFiles}  DESTINATION include/${target_name})

solvers/grainGranFlow/createDEMComponents.hpp

@@ -17,10 +17,17 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
+REPORT(0)<<"Reading shape dictionary ..."<<END_REPORT;
+pFlow::grainShape grains
+(
+	pFlow::shapeFile__,
+	&Control.caseSetup(),
+	proprties
+);
 
 //
 REPORT(0)<<"\nReading sphere particles . . ."<<END_REPORT;
-pFlow::grainParticles grnParticles(Control, proprties);
+pFlow::grainParticles grnParticles(Control, grains);
 
 //
 REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

solvers/iterateSphereParticles/createDEMComponents.hpp

@@ -18,8 +18,16 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+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, proprties);
+pFlow::sphereParticles sphParticles(Control, spheres);
 
 WARNING<<"Particle insertion has not been set yet!"<<END_WARNING;

solvers/sphereGranFlow/createDEMComponents.hpp

@@ -17,10 +17,18 @@ 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, proprties);
+
+
+pFlow::sphereParticles sphParticles(Control, spheres);
 
 //
 REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;

src/Geometry/geometryMotion/geometryMotions.cpp

@@ -26,4 +26,6 @@ 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

@@ -24,6 +24,7 @@ Licence:
 #include "geometryMotion.hpp"
 #include "stationaryWall.hpp"
 #include "rotatingAxisMotion.hpp"
+#include "conveyorBeltMotion.hpp"
 //#include "multiRotatingAxisMotion.hpp"
 #include "vibratingMotion.hpp"
 
@@ -37,6 +38,8 @@ using rotationAxisMotionGeometry = geometryMotion<rotatingAxisMotion>;
 
 using stationaryGeometry = geometryMotion<stationaryWall>;
 
+using conveyorBeltMotionGeometry = geometryMotion<conveyorBeltMotion>;
+
 //typedef geometryMotion<multiRotatingAxisMotion> multiRotationAxisMotionGeometry;
 
 

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

@@ -37,7 +37,8 @@ bool intAllActive(
 	real dt, 
 	realx3Field_D& y, 
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1)
+	realx3PointField_D& dy1,
+	real damping = 1.0)
 {
 
 	auto d_dy = dy.deviceView();
@@ -49,7 +50,7 @@ bool intAllActive(
 		"AdamsBashforth2::correct",
 		rpIntegration (activeRng.start(), activeRng.end()),
 		LAMBDA_HD(uint32 i){
-			d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+			d_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();
@@ -62,7 +63,8 @@ bool intScattered
 	real dt, 
 	realx3Field_D& y,
 	realx3PointField_D& dy,
-	realx3PointField_D& dy1
+	realx3PointField_D& dy1,
+	real damping = 1.0
 )
 {
 
@@ -78,7 +80,7 @@ bool intScattered
 		LAMBDA_HD(uint32 i){
 			if( activeP(i))
 			{
-				d_y[i] +=  dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
+				d_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];
 			}
 		});
@@ -112,8 +114,11 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 		zero3,
 		zero3
 	),
+	initialValField_(initialValField),
 	boundaryList_(pStruct, method, *this)
-{}
+{
+	realx3PointField_D::addEvent(message::ITEMS_INSERT);
+}
 
 void pFlow::AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested()
 {
@@ -142,18 +147,19 @@ bool pFlow::AdamsBashforth2::correct
 (
 	real dt,
 	realx3PointField_D& y,
-	realx3PointField_D& dy
+	realx3PointField_D& dy,
+	real damping
 )
 {
 	auto& dy1l = dy1();
 	bool success = false;
 	if(dy1l.isAllActive())
 	{
-		success = intAllActive(dt, y.field(), dy, dy1l);
+		success = intAllActive(dt, y.field(), dy, dy1(), damping);
 	}
 	else
 	{
-		success = intScattered(dt, y.field(), dy, dy1l);
+		success = intScattered(dt, y.field(), dy, dy1(), damping);
 	}
 
 	success = success && boundaryList_.correct(dt, y, dy);
@@ -171,11 +177,11 @@ bool pFlow::AdamsBashforth2::correctPStruct(
 	bool success = false;
 	if(dy1l.isAllActive())
 	{
-		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1l);
+		success = intAllActive(dt, pStruct.pointPosition(), vel,  dy1());
 	}
 	else
 	{
-		success = intScattered(dt, pStruct.pointPosition(), vel, dy1l);
+		success = intScattered(dt, pStruct.pointPosition(), vel,  dy1());
 	}
 
 	success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
@@ -183,11 +189,21 @@ bool pFlow::AdamsBashforth2::correctPStruct(
 	return success;
 }
 
-
+/*bool pFlow::AdamsBashforth2::hearChanges
-bool pFlow::AdamsBashforth2::setInitialVals(
+(
-	const int32IndexContainer& newIndices,
+	const timeInfo &ti, 
-	const realx3Vector& y)
+	const message &msg, 
+	const anyList &varList
+)
 {
-	return true;
+	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

@@ -41,10 +41,14 @@ class AdamsBashforth2
 {
 private:
 
+	const realx3Field_D& initialValField_;
+
 	boundaryIntegrationList boundaryList_;
 
 	friend class processorAB2BoundaryIntegration;
 
+protected:
+	
 	const auto& dy1()const
 	{
 		return static_cast<const realx3PointField_D&>(*this);
@@ -55,6 +59,16 @@ private:
 		return static_cast<realx3PointField_D&>(*this);
 	}
 
+	auto& initialValField()
+	{
+		return initialValField_;
+	}
+
+	boundaryIntegrationList& boundaryList()
+	{
+		return boundaryList_;
+	}
+	
 public:
 
 	/// Class info
@@ -70,7 +84,7 @@ public:
 			const realx3Field_D& initialValField);
 
 		/// Destructor 
-		~AdamsBashforth2()final = default;
+		~AdamsBashforth2()override = default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -102,32 +116,21 @@ public:
 		bool correct(
 			real dt, 
 			realx3PointField_D& y, 
-			realx3PointField_D& dy) final;
+			realx3PointField_D& dy,
+			real damping = 1.0) override;
 
 		bool correctPStruct(
 			real dt, 
 			pointStructure& pStruct, 
-			realx3PointField_D& vel) final;
+			realx3PointField_D& vel) override;
-			
 		
 		/*bool hearChanges
 		(
-			real t,
+			const timeInfo& ti,
-			real dt,
-			uint32 iter,
 			const message& msg, 
 			const anyList& varList
 		) override;*/
 
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) final;
-
-		bool needSetInitialVals()const final
-		{
-			return false;
-		}
-
 };
 
 

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -20,90 +20,184 @@ Licence:
 
 #include "AdamsBashforth3.hpp"
 
-//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
+namespace pFlow
-
-pFlow::AdamsBashforth3::AdamsBashforth3
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB3History>>(
-			objectFile(
-				groupNames(baseName,"AB3History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB3History({zero3,zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth3::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto activeRng = dy1.activeRange();
-	{
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth3::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth3::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){
+		LAMBDA_HD(uint32 i){
-			auto ldy = d_dy[i];
+			d_y[i] += damping* dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
-			d_y[i] += dt*( static_cast<real>(23.0 / 12.0) * ldy 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
-				- static_cast<real>(16.0 / 12.0) * d_history[i].dy1_ 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]) ;
-				+ static_cast<real>(5.0 / 12.0) * d_history[i].dy2_);
+			
-			d_history[i] = {ldy ,d_history[i].dy1_};
+			d_dy2[i] = 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,
+	realx3PointField_D& dy2,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth3::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping * dt*( static_cast<real>(23.0 / 12.0) * d_dy[i] 
+						- static_cast<real>(16.0 / 12.0) * d_dy1[i] 
+						+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
+			
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
+
+pFlow::AdamsBashforth3::AdamsBashforth3
+(
+	const word& baseName,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField
+)
+:
+	AdamsBashforth2(baseName, pStruct, method, initialValField),
+	dy2_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy2"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+
+{
+	
+}
+
+void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
+	dy2_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+
+
+bool pFlow::AdamsBashforth3::correct
+(
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth3::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth3::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -22,48 +22,14 @@ Licence:
 #define __AdamsBashforth3_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth2.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB3History
-{
-	TypeInfoNV("AB3History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB3History& ab3)
-{
-	str.readBegin("AB3History");
-
-    str >> ab3.dy1_;
-    str >> ab3.dy2_;
-    
-    str.readEnd("AB3History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB3History& ab3)
-{
-	str << token::BEGIN_LIST << ab3.dy1_
-	    << token::SPACE << ab3.dy2_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Third order Adams-Bashforth integration method for solving ODE
@@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
  */
 class AdamsBashforth3
 :
-	public integration
+	public AdamsBashforth2
 {
+private:
+	
+	friend class processorAB3BoundaryIntegration;
+
+	realx3PointField_D 		dy2_;
+
 protected:
 
-	/// Integration history
+	const auto& dy2()const
-	pointField<VectorSingle,AB3History>& history_;
+	{
+		return dy2_;
+	}
+
+	auto& dy2()
+	{
+		return dy2_;
+	}
 	
-	/// Range policy for integration kernel
-	using rpIntegration = Kokkos::RangePolicy<
-			DefaultExecutionSpace,
-			Kokkos::Schedule<Kokkos::Static>,
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -96,17 +69,13 @@ public:
 		/// Construct from components
 		AdamsBashforth3(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth3>(*this);
-		}
 
 		/// Destructor 
-		virtual ~AdamsBashforth3()=default;
+		~AdamsBashforth3() override =default;
 
 		/// Add this to the virtual constructor table 
 		add_vCtor(
@@ -117,43 +86,39 @@ public:
 
 	// - Methods
 
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
-			realx3Vector_D & UNUSED(y),
-			realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, 
+		/// return integration method 
-			realx3Vector_D & y,
+		word method()const override
-			realx3Vector_D& dy) override;
-
-		bool setInitialVals(
-			const int32IndexContainer& newIndices,
-			const realx3Vector& y) override;
-
-		bool needSetInitialVals()const override
 		{
-			return false;
+			return "AdamsBashforth3";
 		}
 
-		/// Integrate on all points in the active range 
-		bool intAll(
-			real dt, 
-			realx3Vector_D& y, 
-			realx3Vector_D& dy, 
-			range activeRng);
 		
-		/// Integrate on active points in the active range 
+
-		template<typename activeFunctor>
+		bool correct(
-		bool intRange(
 			real dt, 
-			realx3Vector_D& y,
+			realx3PointField_D& y, 
-			realx3Vector_D& dy,
+			realx3PointField_D& dy,
-			activeFunctor activeP );
+			real damping = 1.0) override;
+
+		bool correctPStruct(
+			real dt, 
+			pointStructure& pStruct, 
+			realx3PointField_D& vel) override;
+			
+		
+		/*bool hearChanges
+		(
+			const timeInfo& ti,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth3::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -181,7 +146,7 @@ bool pFlow::AdamsBashforth3::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -20,96 +20,186 @@ Licence:
 
 #include "AdamsBashforth4.hpp"
 
-
+namespace pFlow
-
-pFlow::AdamsBashforth4::AdamsBashforth4
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB4History>>(
-			objectFile(
-				groupNames(baseName,"AB4History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB4History({zero3,zero3, zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth4::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto d_dy3= dy3.deviceView();
-	{
+	auto activeRng = dy1.activeRange();
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth4::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){				
+		LAMBDA_HD(uint32 i){
-			d_y[i] += dt*( 
+			d_y[i] += damping * dt*( 
-						static_cast<real>(55.0 / 24.0) * d_dy[i]
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
-					- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
-					+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
-					- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i]);
-					);
-			d_history[i].dy3_ = d_history[i].dy2_;
-			d_history[i].dy2_ = d_history[i].dy1_;
-			d_history[i].dy1_ = d_dy[i];
-			
 
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = 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,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth4::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(55.0 / 24.0) * d_dy[i]
+				- static_cast<real>(59.0 / 24.0) * d_dy1[i] 
+				+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
+				- static_cast<real>( 9.0 / 24.0) * d_dy3[i] );
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth4::AdamsBashforth4
+(
+	const word& baseName,
+	pointStructure& pStruct,
+	const word& method,
+	const realx3Field_D& initialValField
+)
+:
+	AdamsBashforth3(baseName, pStruct, method, initialValField),
+	dy3_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy3"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
+	dy3_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth4::correct
+(
+	real dt, 
+	realx3PointField_D& y, 
+	realx3PointField_D& dy,
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth4::correctPStruct(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel)
+{
+	
+	bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+/*bool pFlow::AdamsBashforth4::hearChanges
+(
+	const timeInfo &ti, 
+	const message &msg, 
+	const anyList &varList
+)
+{
+	if(msg.equivalentTo(message::ITEMS_INSERT))
+	{
+		return insertValues(varList, initialValField().deviceViewAll(), dy1()) &&
+		insertValues(varList, initialValField().deviceViewAll(), dy2()) && 
+		insertValues(varList, initialValField().deviceViewAll(), dy3());
+	}
+	else
+	{
+		return realx3PointField_D::hearChanges(ti, msg, varList);
+	}
+}*/

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -22,53 +22,14 @@ Licence:
 #define __AdamsBashforth4_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth3.hpp"
 #include "pointFields.hpp"
+#include "boundaryIntegrationList.hpp"
+
 
 namespace pFlow
 {
 
-struct AB4History
-{
-	TypeInfoNV("AB4History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB4History& ab4)
-{
-	str.readBegin("AB4History");
-
-    str >> ab4.dy1_;
-    str >> ab4.dy2_;
-    str >> ab4.dy3_;
-    
-    str.readEnd("AB4History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB4History& ab4)
-{
-	str << token::BEGIN_LIST << ab4.dy1_
-	    << token::SPACE << ab4.dy2_
-	    << token::SPACE << ab4.dy3_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
-
 /**
  * Fourth order Adams-Bashforth integration method for solving ODE
  * 
@@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
  */
 class AdamsBashforth4
 :
-	public integration
+	public AdamsBashforth3
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy3_;
+
 protected:
 
-	/// Integration history
+	const auto& dy3()const
-	pointField<VectorSingle,AB4History>& history_;
+	{
+		return dy3_;
+	}
 
-	/// Range policy for integration kernel
+	auto& dy3()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return dy3_;
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -100,17 +67,14 @@ public:
 		/// Construct from components
 		AdamsBashforth4(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField);
+		
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth4>(*this);
-		}
 
 		/// Destructor
-		virtual ~AdamsBashforth4()=default;
+		~AdamsBashforth4() override =default;
 
 		/// Add a this to the virtual constructor table 
 		add_vCtor(
@@ -121,77 +85,36 @@ public:
 
 	// - Methods
 		
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
+
-			realx3Vector_D & UNUSED(y),
+		/// return integration method 
-			realx3Vector_D& UNUSED(dy)) override;
+		word method()const override
+		{
+			return "AdamsBashforth4";
+		}
 
 		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::AdamsBashforth4::intRange(
-	real dt, 
-	realx3Vector_D& y,
-	realx3Vector_D& dy,
-	activeFunctor activeP )
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
-	auto activeRng = activeP.activeRange();
 
-	Kokkos::parallel_for(
-		"AdamsBashforth4::correct",
-		rpIntegration (activeRng.first, activeRng.second),
-		LAMBDA_HD(int32 i){
-			if( activeP(i))
-			{
-				
-				d_y[i] += dt*( 
-							static_cast<real>(55.0 / 24.0) * d_dy[i]
-						- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_ 
-						+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
-						- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
-						);
-				d_history[i].dy3_ = d_history[i].dy2_;
-				d_history[i].dy2_ = d_history[i].dy1_;
-				d_history[i].dy1_ = d_dy[i];
-			}
-		});
-	Kokkos::fence();
-
-	return true;
-}
 
 } // pFlow
 

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -20,93 +20,178 @@ Licence:
 
 #include "AdamsBashforth5.hpp"
 
-
+namespace pFlow
-
-pFlow::AdamsBashforth5::AdamsBashforth5
-(
-	const word& baseName,
-	repository& owner,
-	const pointStructure& pStruct,
-	const word& method
-)
-:
-	integration(baseName, owner, pStruct, method),
-	history_(
-		owner.emplaceObject<pointField<VectorSingle,AB5History>>(
-			objectFile(
-				groupNames(baseName,"AB5History"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			AB5History({zero3,zero3, zero3})))
-
 {
 
-}
+/// Range policy for integration kernel (alias)
+using rpIntegration = Kokkos::RangePolicy<
+		DefaultExecutionSpace,
+		Kokkos::Schedule<Kokkos::Static>,
+		Kokkos::IndexType<uint32>
+		>;
 
-bool pFlow::AdamsBashforth5::predict
+bool intAllActive(
-(
-	real UNUSED(dt),
-	realx3Vector_D& UNUSED(y),
-	realx3Vector_D& UNUSED(dy)
-)
-{
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::correct
-(
 	real dt, 
-	realx3Vector_D& y,
+	realx3Field_D& y, 
-	realx3Vector_D& dy
+	realx3PointField_D& dy,
-)
+	realx3PointField_D& dy1,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0)
 {
 
-	if(this->pStruct().allActive())
+	auto d_dy = dy.deviceView();
-	{
+	auto d_y  = y.deviceView();
-		return intAll(dt, y, dy, this->pStruct().activeRange());
+	auto d_dy1= dy1.deviceView();
-	}
+	auto d_dy2= dy2.deviceView();
-	else
+	auto d_dy3= dy3.deviceView();
-	{
+	auto d_dy4= dy4.deviceView();
-		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+	auto activeRng = dy1.activeRange();
-	}
-
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::setInitialVals(
-	const int32IndexContainer& newIndices,
-	const realx3Vector& y)
-{
-	return true;
-}
-
-bool pFlow::AdamsBashforth5::intAll(
-	real dt, 
-	realx3Vector_D& y, 
-	realx3Vector_D& dy, 
-	range activeRng)
-{
-	auto d_dy = dy.deviceViewAll();
-	auto d_y  = y.deviceViewAll();
-	auto d_history = history_.deviceViewAll();
 
 	Kokkos::parallel_for(
 		"AdamsBashforth5::correct",
-		rpIntegration (activeRng.first, activeRng.second),
+		rpIntegration (activeRng.start(), activeRng.end()),
-		LAMBDA_HD(int32 i){				
+		LAMBDA_HD(uint32 i){
-			d_y[i] += dt*( 
+			d_y[i] += damping*  dt*( 
 				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
-				- static_cast<real>(2774.0 / 720.0) * d_history[i].dy1_ 
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
-				+ static_cast<real>(2616.0 / 720.0) * d_history[i].dy2_
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
-				- static_cast<real>(1274.0 / 720.0) * d_history[i].dy3_
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
-				+ static_cast<real>( 251.0 / 720.0) * d_history[i].dy4_
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
-				);
+				
-			d_history[i] = {d_dy[i] ,d_history[i].dy1_, d_history[i].dy2_, d_history[i].dy3_};
+			d_dy4[i] = d_dy3[i];
+			d_dy3[i] = d_dy2[i];
+			d_dy2[i] = 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,
+	realx3PointField_D& dy2,
+	realx3PointField_D& dy3,
+	realx3PointField_D& dy4,
+	real damping = 1.0
+)
+{
+
+	auto d_dy 		= dy.deviceView();
+	auto d_y  		= y.deviceView();
+	auto d_dy1		= dy1.deviceView();
+	auto d_dy2		= dy2.deviceView();
+	auto d_dy3		= dy3.deviceView();
+	auto d_dy4		= dy4.deviceView();
+	auto activeRng 	= dy1.activeRange();
+	const auto& activeP = dy1.activePointsMaskDevice();
+
+	Kokkos::parallel_for(
+		"AdamsBashforth5::correct",
+		rpIntegration (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			if( activeP(i))
+			{
+				d_y[i] += damping*  dt*( 
+				  static_cast<real>(1901.0 / 720.0) * d_dy[i]
+				- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+				+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
+				- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+				+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
+
+				d_dy4[i] = d_dy3[i];
+				d_dy3[i] = d_dy2[i];
+				d_dy2[i] = d_dy1[i];
+				d_dy1[i] = d_dy[i];
+			}
+		});
+	Kokkos::fence();
+
+	return true;
+}
+
+}
+
+pFlow::AdamsBashforth5::AdamsBashforth5
+(
+	const word &baseName, 
+	pointStructure &pStruct, 
+	const word &method, 
+	const realx3Field_D &initialValField
+)
+:
+	AdamsBashforth4(baseName, pStruct, method, initialValField),
+	dy4_
+	(
+		objectFile
+		(
+			groupNames(baseName,"dy4"),
+			pStruct.time().integrationFolder(),
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		pStruct,
+		zero3,
+		zero3
+	)
+{}
+
+void pFlow::AdamsBashforth5::updateBoundariesSlaveToMasterIfRequested()
+{
+	AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested();
+	dy4_.updateBoundariesSlaveToMasterIfRequested();
+}
+
+bool pFlow::AdamsBashforth5::correct
+(
+	real dt, 
+	realx3PointField_D &y, 
+	realx3PointField_D &dy, 
+	real damping
+)
+{
+	
+	bool success = false;
+	if(y.isAllActive())
+	{
+		success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+	else
+	{
+		success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
+	}
+
+	success = success && boundaryList().correct(dt, y, dy);
+
+	return success;
+}
+
+bool pFlow::AdamsBashforth5::correctPStruct
+(
+	real dt, 
+	pointStructure &pStruct, 
+	realx3PointField_D &vel
+)
+{
+bool success = false;
+	if(dy2().isAllActive())
+	{
+		success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+	else
+	{
+		success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
+	}
+
+	success = success && boundaryList().correctPStruct(dt, pStruct, vel);
+
+	return success;
+}
+
+

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -22,54 +22,12 @@ Licence:
 #define __AdamsBashforth5_hpp__
 
 
-#include "integration.hpp"
+#include "AdamsBashforth4.hpp"
 #include "pointFields.hpp"
 
 namespace pFlow
 {
 
-struct AB5History
-{
-	TypeInfoNV("AB5History");
-
-	realx3 dy1_={0,0,0};
-	realx3 dy2_={0,0,0};
-	realx3 dy3_={0,0,0};
-	realx3 dy4_={0,0,0};
-};
-
-
-INLINE_FUNCTION
-iIstream& operator>>(iIstream& str, AB5History& ab5)
-{
-	str.readBegin("AB5History");
-
-    str >> ab5.dy1_;
-    str >> ab5.dy2_;
-    str >> ab5.dy3_;
-    str >> ab5.dy4_;
-    
-    str.readEnd("AB5History");
-
-    str.check(FUNCTION_NAME);
-
-	return str;
-	
-}
-
-INLINE_FUNCTION
-iOstream& operator<<(iOstream& str, const AB5History& ab5)
-{
-	str << token::BEGIN_LIST << ab5.dy1_
-	    << token::SPACE << ab5.dy2_
-	    << token::SPACE << ab5.dy3_
-	    << token::SPACE << ab5.dy4_
-		<< token::END_LIST;	    
-
-	str.check(FUNCTION_NAME);
-
-	return str;
-}
 
 /**
  * Fifth order Adams-Bashforth integration method for solving ODE
@@ -78,19 +36,25 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
  */
 class AdamsBashforth5
 :
-	public integration
+	public AdamsBashforth4
 {
+private:
+
+	friend class processorAB4BoundaryIntegration;
+
+	realx3PointField_D 		dy4_;
+
 protected:
 
-	/// Integration history
+	const auto& dy4()const
-	pointField<VectorSingle,AB5History>& history_;
+	{
+		return dy4_;
+	}
 
-	/// Range policy for integration kernel
+	auto& dy4()
-	using rpIntegration = Kokkos::RangePolicy<
+	{
-			DefaultExecutionSpace,
+		return dy4_;
-			Kokkos::Schedule<Kokkos::Static>,
+	}
-			Kokkos::IndexType<int32>
-			>;
 
 public:
 
@@ -99,20 +63,19 @@ public:
 
 	// - Constructors
 		
+		/// Construct from components
 		AdamsBashforth5(
 			const word& baseName,
-			repository& owner,
+			pointStructure& pStruct,
-			const pointStructure& pStruct,
+			const word& method,
-			const word& method);
+			const realx3Field_D& initialValField);
 		
-		uniquePtr<integration> clone()const override
-		{
-			return makeUnique<AdamsBashforth5>(*this);
-		}
 		
-		virtual ~AdamsBashforth5()=default;
 
-		/// Add this to the virtual constructor table 
+		/// Destructor
+		~AdamsBashforth5() override =default;
+
+		/// Add a this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth5,
@@ -121,44 +84,29 @@ public:
 
 	// - Methods
 		
-		bool predict(
+		void updateBoundariesSlaveToMasterIfRequested()override;
-			real UNUSED(dt),
+
-			realx3Vector_D & UNUSED(y),
+		/// return integration method 
-			realx3Vector_D& UNUSED(dy)) override;
+		word method()const override
+		{
+			return "AdamsBashforth5";
+		}
 
 		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 
-		template<typename activeFunctor>
-		bool intRange(
-			real dt,
-			realx3Vector_D& y,
-			realx3Vector_D& dy,
-			activeFunctor activeP );
 
 };
 
 
-template<typename activeFunctor>
+/*template<typename activeFunctor>
 bool pFlow::AdamsBashforth5::intRange(
 	real dt, 
 	realx3Vector_D& y,
@@ -190,7 +138,7 @@ bool pFlow::AdamsBashforth5::intRange(
 	Kokkos::fence();
 
 	return true;
-}
+}*/
 
 } // pFlow
 

src/Integration/CMakeLists.txt

@@ -4,9 +4,10 @@ integration/integration.cpp
 boundaries/boundaryIntegration.cpp
 boundaries/boundaryIntegrationList.cpp 
 AdamsBashforth2/AdamsBashforth2.cpp
-#AdamsBashforth5/AdamsBashforth5.cpp
+AdamsBashforth3/AdamsBashforth3.cpp
-#AdamsBashforth4/AdamsBashforth4.cpp
+AdamsBashforth4/AdamsBashforth4.cpp
-#AdamsBashforth3/AdamsBashforth3.cpp
+AdamsBashforth5/AdamsBashforth5.cpp
+
 #AdamsMoulton3/AdamsMoulton3.cpp
 #AdamsMoulton4/AdamsMoulton4.cpp
 #AdamsMoulton5/AdamsMoulton5.cpp

src/Integration/boundaries/boundaryIntegration.hpp

@@ -51,9 +51,7 @@ public:
     );
 
     bool hearChanges(
-        real t,
+        const timeInfo& ti,
-        real dt,
-        uint32 iter,
         const message &msg,
         const anyList &varList) override
     {
@@ -77,6 +75,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+    
     static
     uniquePtr<boundaryIntegration> create(
         const boundaryBase& boundary,

src/Integration/boundaries/boundaryIntegrationList.cpp

@@ -7,11 +7,11 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
     integration &intgrtn
 )
 :
-    ListPtr<boundaryIntegration>(6),
+    boundaryListPtr<boundaryIntegration>(),
 	boundaries_(pStruct.boundaries())
 {
 	
-    for(uint32 i=0; i<6; i++)
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,
@@ -19,25 +19,23 @@ pFlow::boundaryIntegrationList::boundaryIntegrationList(
 				boundaries_[i],
 				pStruct, 
 				method,
-                intgrtn
+				intgrtn
-            )
+			)
-        );
+		);
 	}	
-
 }
 
 bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
 {
-	for(auto& bndry:*this)
+	ForAllActiveBoundariesPtr(i,this)
 	{
-		if(!bndry->correct(dt, y, dy))
+		if(!boundaryPtr(i)->correct(dt, y, dy))
 		{
 			fatalErrorInFunction<<"Error in correcting boundary "<<
-			bndry->boundaryName()<<endl;
+			boundaryPtr(i)->boundaryName()<<endl;
 			return false;
 		}
-		
+	}
-    }
     return true;
 }
 
@@ -46,14 +44,15 @@ bool pFlow::boundaryIntegrationList::correctPStruct(
 	pointStructure &pStruct, 
 	const realx3PointField_D &vel)
 {
-	for(auto& bndry:*this)
+	ForAllActiveBoundariesPtr(i,this)
 	{
-		if(!bndry->correctPStruct(dt, vel))
+		if(!boundaryPtr(i)->correctPStruct(dt, vel))
 		{
 			fatalErrorInFunction<<"Error in correcting boundary "<<
-			bndry->boundaryName()<<" in pointStructure."<<endl;
+			boundaryPtr(i)->boundaryName()<<" in pointStructure."<<endl;
 			return false;
 		}
 	}
+	
     return true;
 }

src/Integration/boundaries/boundaryIntegrationList.hpp

@@ -4,7 +4,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryIntegration.hpp"
 
 
@@ -15,7 +15,7 @@ class integration;
 
 class boundaryIntegrationList
 :
-    public ListPtr<boundaryIntegration>
+    public boundaryListPtr<boundaryIntegration>
 {
 private:
 

src/Integration/integration/integration.cpp

@@ -22,17 +22,39 @@ Licence:
 #include "pointStructure.hpp"
 #include "repository.hpp"
 
-pFlow::integration::integration
+bool pFlow::integration::insertValues
 (
-	const word& baseName,
+	const anyList& varList, 
-	pointStructure& pStruct,
+	const deviceViewType1D<realx3>& srcVals,
-	const word&,
+	realx3PointField_D& dstFeild
-	const realx3Field_D&
 )
-:
+{
-	owner_(*pStruct.owner()),
+	const word eventName = message::eventName(message::ITEMS_INSERT);
-	pStruct_(pStruct),
+
-	baseName_(baseName)	
+	if(!varList.checkObjectType<uint32IndexContainer>(eventName))
+	{
+		fatalErrorInFunction<<"Insertion failed in integration for "<< baseName_<<
+		", variable "<< eventName <<
+		" does not exist or the type is incorrect"<<endl;
+		return false;
+	}
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		eventName);
+	
+	dstFeild.field().insertSetElement(indices, srcVals);
+
+    return true;
+}
+
+pFlow::integration::integration(
+    const word &baseName,
+    pointStructure &pStruct,
+    const word &,
+    const realx3Field_D &)
+    : owner_(*pStruct.owner()),
+      pStruct_(pStruct),
+      baseName_(baseName)
 {}
 
 

src/Integration/integration/integration.hpp

@@ -63,6 +63,13 @@ private:
 		/// The base name for integration 
 		const word baseName_;
 
+protected:
+
+	bool insertValues(
+		const anyList& varList, 
+		const deviceViewType1D<realx3>& srcVals,
+		realx3PointField_D& dstFeild);
+
 public:
 
 	/// Type info
@@ -146,22 +153,11 @@ public:
 
 		/// Correction/main integration step
 		virtual 
-		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
+		bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
 
 		virtual 
 		bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;
 
-		/// Set the initial values for new indices 
-		virtual 
-		bool setInitialVals(
-				const int32IndexContainer& newIndices,
-				const realx3Vector& y) = 0;
-
-		/// Check if the method requires any set initial vals
-		virtual 
-		bool needSetInitialVals()const = 0;
-
-
 	/// Create the polymorphic object based on inputs
 	static
 		uniquePtr<integration> create(

src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp

@@ -25,12 +25,6 @@ Licence:
 #include "symArrays.hpp"
 
 
-
-
-
-
-
-
 namespace pFlow::cfModels
 {
 
@@ -159,7 +153,7 @@ protected:
 		{
 
 			etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
-					sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
+					sqrt(log(pow(et[i],2))+ pow(Pi,2));
 		}
 
 		Vector<linearProperties> prop("prop", nElem);
@@ -285,8 +279,8 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 
@@ -299,29 +293,24 @@ public:
 		else if (addDissipationModel_==3)
 		{
 			auto pie =3.14;
-			prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
+			prop.en_ = exp((pow(f_,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
 		}
 
-
-
-
 		real ethan_  = -2.0*log(prop.en_)*sqrt(prop.kn_)/
-					sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(prop.en_),2)+ pow(Pi,2));
 
 			//REPORT(0)<<"\n en n is  :  "<<END_REPORT;
 			//REPORT(0)<< prop.en_ <<END_REPORT;
 
-		FCn = (  -pow(f_,3.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,1.5) * ethan_ * vrn)*Nij;
+		FCn = (  -pow(f_,3)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,static_cast<real>(1.5)) * ethan_ * vrn)*Nij;
-		FCt = ( -pow(f_,3.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,1.5) *  prop.ethat_*Vt);
+		FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast<real>(1.5)) *  prop.ethat_*Vt);
-		
-
 				
 		real ft = length(FCt);
 		real ft_fric = prop.mu_ * length(FCn);
 		
 		if(ft > ft_fric)
 		{
-			if( length(history.overlap_t_) >static_cast<real>(0.0))
+			if( length(history.overlap_t_) >zero)
 			{
 				if constexpr (limited)
 				{

src/Interaction/Models/contactForce/cGNonLinearCF.hpp

@@ -0,0 +1,307 @@
+/*------------------------------- 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 __cGNonLinearCF_hpp__
+#define __cGNonLinearCF_hpp__
+
+#include "types.hpp"
+
+namespace pFlow::cfModels
+{
+
+template<bool limited=true>
+class cGNonLinear
+{
+public:
+
+	struct contactForceStorage
+	{
+		realx3 overlap_t_ = 0.0;
+	};
+
+	struct cGNonLinearProperties
+	{		
+		real 	Yeff_	= 1000000.0;
+		real 	Geff_	= 8000000.0;
+		real 	en_ 	= 1.0;
+		real  	mu_ 	= 0.00001;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(){}
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
+			Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
+		{}		
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
+
+		INLINE_FUNCTION_HD
+		~cGNonLinearProperties() = default;
+	};
+
+protected:
+
+	using cGNonLinearArrayType = 	symArray<cGNonLinearProperties>;
+
+	int32 						numMaterial_ = 0;
+
+	ViewType1D<real>  		 	rho_;
+
+	int32 						addDissipationModel_;
+
+	cGNonLinearArrayType  		nonlinearProperties_;
+
+	bool readNonLinearDictionary(const dictionary& dict)
+	{
+		auto Yeff = dict.getVal<realVector>("Yeff");
+		auto Geff = dict.getVal<realVector>("Geff");
+		auto nu = dict.getVal<realVector>("nu");
+		auto en = dict.getVal<realVector>("en");
+		auto mu = dict.getVal<realVector>("mu");
+
+		auto nElem = Yeff.size();
+
+		if(nElem != nu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
+			return false;
+		}
+
+		if(nElem != en.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
+			return false;
+		}
+
+
+		if(nElem != mu.size())
+		{
+			fatalErrorInFunction<<
+			"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
+			return false;
+		}
+
+		
+		// check if size of vector matchs a symetric array
+		uint32 nMat;
+		if( !cGNonLinearArrayType::getN(nElem, nMat) )
+		{
+			fatalErrorInFunction<<
+			"sizes of properties do not match a symetric array with size ("<<
+			numMaterial_<<"x"<<numMaterial_<<").\n";
+			return false;
+		}
+		else if( numMaterial_ != nMat)
+		{
+			fatalErrorInFunction<<
+			"size mismatch for porperties. \n";
+			return false;
+		}
+
+		Vector<cGNonLinearProperties> prop("prop",nElem);
+		ForAll(i,Yeff)
+		{
+			prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
+		}
+
+		nonlinearProperties_.assign(prop);
+
+		auto adm = dict.getVal<word>("additionalDissipationModel");
+		if(adm == "none")
+		{
+			addDissipationModel_ = 1;
+		}
+		else if(adm == "LU")
+		{
+			addDissipationModel_ = 2;
+		}
+		else if (adm == "GB")
+		{
+			addDissipationModel_ = 3;
+		}
+		else
+		{
+			addDissipationModel_ = 1;
+		}
+
+		return true;
+
+	}
+
+	static const char* modelName()
+	{
+		if constexpr (limited)
+		{
+			return "cGNonLinearLimited";
+		}
+		else
+		{
+			return "cGNonLinearNonLimited";
+		}
+		return "";
+	}
+
+public:
+
+	TypeInfoNV(modelName());
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(){}
+
+	cGNonLinear(
+		int32 nMaterial,
+		const ViewType1D<real>& rho,
+		const dictionary& dict)
+	:
+		numMaterial_(nMaterial),
+		rho_("rho",nMaterial),
+		nonlinearProperties_("cGNonLinearProperties",nMaterial)
+	{
+
+		Kokkos::deep_copy(rho_,rho);
+		if(!readNonLinearDictionary(dict))
+		{
+			fatalExit;
+		}
+	}
+	
+	INLINE_FUNCTION_HD
+	cGNonLinear(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear(cGNonLinear&&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(const cGNonLinear&) = default;
+
+	INLINE_FUNCTION_HD
+	cGNonLinear& operator=(cGNonLinear&&) = default;
+
+
+	INLINE_FUNCTION_HD
+	~cGNonLinear()=default;
+
+	INLINE_FUNCTION_HD
+	int32 numMaterial()const
+	{
+		return numMaterial_;
+	}
+
+	//// - Methods
+
+	INLINE_FUNCTION_HD
+	void contactForce
+	(
+		const real dt,
+		const uint32 i,
+		const uint32 j,
+		const uint32 propId_i,
+		const uint32 propId_j,
+		const real Ri,
+		const real Rj,
+		const real cGFi,
+		const real cGFj,
+		const real ovrlp_n,
+		const realx3& Vr,
+		const realx3& Nij,
+		contactForceStorage& history,
+		realx3& FCn,
+		realx3& FCt
+	)const
+	{
+
+		const auto prop = nonlinearProperties_(propId_i,propId_j);
+
+		const real f = 2/( 1/cGFi + 1/cGFj );
+
+		real vrn = dot(Vr, Nij);	
+		realx3 Vt = Vr - vrn*Nij;
+
+		history.overlap_t_ += Vt*dt;
+
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
+		real Reff = 1/(1/Ri + 1/Rj);
+
+		real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
+		real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
+
+		real en = prop.en_;
+		if (addDissipationModel_==2)
+		{
+			en = sqrt(1+((pow(prop.en_,2)-1)*f));
+		}
+		else if (addDissipationModel_==3)
+		{
+			
+			en = exp((pow(f,static_cast<real>(1.5))*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
+		}
+		
+		real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
+
+		real ethan_  = -2.0*log(en)*sqrt(Kn)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
+
+		FCn = ( - Kn*ovrlp_n - 
+			   sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
+
+		FCt = (f*f)*(- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
+
+		real ft = length(FCt);
+		real ft_fric = prop.mu_ * length(FCn);
+
+		// apply friction 
+		if(ft > ft_fric)
+		{
+			if( length(history.overlap_t_) >0.0)
+			{
+				if constexpr (limited)
+				{
+					real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
+					FCt *= (ft_fric/ft);
+                	history.overlap_t_ = - FCt/(f*f*kt);
+				}
+				else
+				{
+					FCt = (FCt/ft)*ft_fric;	
+				}
+				
+			}
+			else
+			{
+				FCt = 0.0;
+			}
+		}
+
+	}
+   
+
+};
+
+} //pFlow::CFModels
+
+#endif

src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp

@@ -26,11 +26,6 @@ Licence:
 
 
 
-
-
-
-
-
 namespace pFlow::cfModels
 {
 
@@ -49,16 +44,15 @@ public:
 	{		
 		real 	kn_ 	= 1000.0;
 		real 	kt_	= 800.0;
-		real 	en_ 	= 0.0;
+		real 	en_ 	= 1.0;
-		real 	ethat_ 	= 0.0;
+		real  mu_ 	= 0.00001;
-		real    mu_ 	= 0.00001;
 
 		INLINE_FUNCTION_HD
 		linearProperties(){}
 
 		INLINE_FUNCTION_HD
-		linearProperties(real kn, real kt, real en, real etha_t, real mu  ):
+		linearProperties(real kn, real kt, real en, real mu  ):
-			kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu)
+			kn_(kn), kt_(kt), en_(en), mu_(mu)
 		{}		
 
 		INLINE_FUNCTION_HD
@@ -93,7 +87,6 @@ protected:
 		auto kn = dict.getVal<realVector>("kn");
 		auto kt = dict.getVal<realVector>("kt");
 		auto en = dict.getVal<realVector>("en");
-		auto et = dict.getVal<realVector>("et");
 		auto mu = dict.getVal<realVector>("mu");
 		
 
@@ -114,12 +107,6 @@ protected:
 			return false;
 		}
 
-		if(nElem != et.size())
-		{
-			fatalErrorInFunction<<
-			"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
-			return false;
-		}
 
 		if(nElem != mu.size())
 		{
@@ -148,30 +135,16 @@ protected:
 		}
 
 
-		realVector etha_t("etha_t", nElem);
-
-
-		ForAll(i , kn)
-		{
-
-
-			etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
-					sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
-					
-					
-		}
-
 		Vector<linearProperties> prop("prop", nElem);
 		ForAll(i,kn)
 		{
-			prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i]  };
+			prop[i] = {kn[i], kt[i], en[i], mu[i]  };
 		}
 
 		linearProperties_.assign(prop);
 		
-				auto adm = dict.getVal<word>("additionalDissipationModel");
-
 
+		auto adm = dict.getVal<word>("additionalDissipationModel");
 		if(adm == "none")
 		{
 			addDissipationModel_ = 1;
@@ -282,27 +255,40 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 		
-		
+		real en = prop.en_;
-				if (addDissipationModel_==2)
+		if (addDissipationModel_==2)
 		{
-			prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_));
+			en = sqrt(1+((pow(prop.en_,2)-1)*f_));
 		}
-				else if (addDissipationModel_==3)
+		else if (addDissipationModel_==3)
 		{
-			auto pie =3.14;
+			
-			prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) ));
+			en = exp(
+				(
+					pow(f_,static_cast<real>(1.5))*
+					log(prop.en_)*
+					sqrt( 
+						(1-((pow(log(prop.en_),static_cast<real>(2.0))
+						)/
+						(
+							pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0)))))/
+							(1-(pow(f_,3)*(pow(log(prop.en_),2))/
+							(pow(log(prop.en_),static_cast<real>(2.0))+
+							pow(Pi,static_cast<real>(2.0))))) ) ));
 		}
-		real ethan_  = -2.0*log(prop.en_)*sqrt(prop.kn_)/
+
-			sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0));
+		real ethan_  = -2.0*log(en)*sqrt(prop.kn_)/
+			sqrt(pow(log(en),2)+ pow(Pi,2));
 		
 
-		FCn = (  -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij;
+		FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij;
-		FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) *  prop.ethat_*Vt);
+		FCt = ( -f_*prop.kt_ * history.overlap_t_);
 		
 
 		

src/Interaction/Models/contactForce/linearCF.hpp

@@ -139,10 +139,10 @@ protected:
 		ForAll(i , kn)
 		{
 			etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/
-					sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(et[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 		}
 
 		Vector<linearProperties> prop("prop", nElem);
@@ -243,8 +243,8 @@ public:
 
 		history.overlap_t_ += Vt*dt;
 
-		real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i];
+		real mi = 3*Pi/4*pow(Ri,static_cast<real>(3.0))*rho_[propId_i];
-		real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j];
+		real mj = 3*Pi/4*pow(Rj,static_cast<real>(3.0))*rho_[propId_j];
 
 		real sqrt_meff = sqrt((mi*mj)/(mi+mj));
 

src/Interaction/Models/contactForce/nonLinearCF.hpp

@@ -131,7 +131,7 @@ protected:
 			// we take out sqrt(meff*K_hertz) here and then consider this term 
 			// when calculating damping part. 
 			etha_n[i] = -2.2664*log(en[i])/
-					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
+					sqrt(pow(log(en[i]),static_cast<real>(2.0))+ pow(Pi,static_cast<real>(2.0)));
 
 			// no damping for tangential part 
 
@@ -255,7 +255,7 @@ public:
 		// apply friction 
 		if(ft > ft_fric)
 		{
-			if( length(history.overlap_t_) >0.0)
+			if( length(history.overlap_t_) >zero)
 			{
 				if constexpr (limited)
 				{

src/Interaction/Models/grainContactForceModels.hpp

@@ -23,6 +23,7 @@ Licence:
 
 #include "cGAbsoluteLinearCF.hpp"
 #include "cGRelativeLinearCF.hpp"
+#include "cGNonLinearCF.hpp"
 
 #include "grainRolling.hpp"
 
@@ -38,6 +39,10 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<fal
 using limitedCGRelativeLinearGrainRolling   = grainRolling<cGRelativeLinear<true>>;
 using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
 
+using limitedCGNonLinearGrainRolling   = grainRolling<cGNonLinear<true>>;
+using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
+
+
 }
 
 

src/Interaction/Models/rolling/grainRolling.hpp

@@ -96,10 +96,10 @@ public:
 		realx3 w_hat = wi-wj;
 		real w_hat_mag = length(w_hat);
 
-		if( !equal(w_hat_mag,0.0) )
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
 			w_hat /= w_hat_mag;
 		else
-			w_hat = 0.0;
+			w_hat = static_cast<real>(0.0);
 
 		auto Reff = (Ri*Rj)/(Ri+Rj);
 

src/Interaction/Models/rolling/normalRolling.hpp

@@ -94,10 +94,10 @@ public:
 		realx3 w_hat = wi-wj;
 		real w_hat_mag = length(w_hat);
 
-		if( !equal(w_hat_mag,0.0) )
+		if( !equal(w_hat_mag,static_cast<real>(0.0)) )
 			w_hat /= w_hat_mag;
 		else
-			w_hat = 0.0;
+			w_hat = static_cast<real>(0.0);
 
 		auto Reff = (Ri*Rj)/(Ri+Rj);
 

src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp

@@ -60,9 +60,9 @@ private:
 		bool force = false) override
 	{
 		
-		const auto& position = Particles().pointPosition().deviceViewAll();
+		auto position = Particles().pointPosition().deviceViewAll();
-		const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
+		auto flags = Particles().dynPointStruct().activePointsMaskDevice();
-		const auto& diam = Particles().boundingSphere().deviceViewAll();
+		auto diam = Particles().boundingSphere().deviceViewAll();
 
 		return ppwContactSearch_().broadSearch(
 			ti.iter(),
@@ -134,7 +134,7 @@ public:
 		ppwContactSearch_ = 
 			makeUnique<SearchMethodType>
 			(
-				dict(),
+				csDict,
 				this->extendedDomainBox(),
 				minD,
 				maxD,
@@ -154,10 +154,8 @@ public:
 		ContactSearch,
 		dictionary);
 
-	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override
+	bool enterBroadSearchBoundary(const timeInfo& ti, bool force=false)const override;
-	{
+	
-		return enterBroadSearch(ti, force) || csBoundaries_.boundariesUpdated();
-	}
 
 	real sizeRatio()const override
 	{
@@ -171,7 +169,12 @@ public:
 	
 };
 
+template <class searchMethod>
+inline bool ContactSearch<searchMethod>::enterBroadSearchBoundary(const timeInfo &ti, bool force) const
+{
+	return performSearch(ti.iter(),force) || csBoundaries_.boundariesUpdated();
 }
 
+}
 
 #endif //__ContactSearch_hpp__

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

@@ -73,9 +73,7 @@ public:
     }
 
     bool hearChanges(
-        real t,
+        const timeInfo& ti,
-        real dt,
-        uint32 iter,
         const message &msg,
         const anyList &varList) override
     {
@@ -83,8 +81,10 @@ public:
         if (msg.equivalentTo(message::BNDR_RESET))
         {
             // do nothing
+            return true;
         }
-        return true;
+        fatalErrorInFunction;
+        return false;
     }
 
     virtual bool broadSearch(
@@ -98,6 +98,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+
     static uniquePtr<boundaryContactSearch> create(
         const dictionary &dict,
         const boundaryBase &boundary,

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.cpp

@@ -1,3 +1,23 @@
+/*------------------------------- 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 "boundaryContactSearchList.hpp"
 #include "boundaryList.hpp"
 
@@ -5,7 +25,7 @@ void pFlow::boundaryContactSearchList::setList(
     const dictionary &dict, 
     const contactSearch &cSearch)
 {
-    for(auto i=0; i<boundaries_.size(); i++)
+    ForAllBoundariesPtr(i, this)
     {
         this->set
         (
@@ -20,7 +40,7 @@ pFlow::boundaryContactSearchList::boundaryContactSearchList(
     const boundaryList& bndrs,
     const contactSearch &cSearch)
     :
-    ListPtr(bndrs.size()), 
+    boundaryListPtr(), 
     boundaries_(bndrs)
 {
     setList(dict, cSearch);

src/Interaction/contactSearch/boundaries/boundaryContactSearchList.hpp

@@ -1,4 +1,24 @@
-#include "ListPtr.hpp"
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "boundaryListPtr.hpp"
 #include "boundaryContactSearch.hpp"
 
 namespace pFlow
@@ -9,7 +29,7 @@ class contactSearch;
 
 class boundaryContactSearchList
 :
-    public ListPtr<boundaryContactSearch>
+    public boundaryListPtr<boundaryContactSearch>
 {
 private:
 

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

@@ -74,6 +74,11 @@ public:
 		csPairContainerType &ppPairs,
 		csPairContainerType &pwPairs,
 		bool force = false) override;
+	
+	bool isActive()const override
+	{
+		return true;
+	}
 };
 
 }

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

@@ -11,7 +11,7 @@ pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch
     const ViewType1D<realx3, memory_space> &normals
 )
 :
-    cellExtent_( max(cellExtent, 0.5 ) ),
+    cellExtent_( max(cellExtent, half ) ),
     numElements_(numElements),
     numPoints_(numPoints),
     vertices_(vertices),

src/Interaction/contactSearch/contactSearch/contactSearch.cpp

@@ -31,11 +31,11 @@ pFlow::contactSearch::contactSearch(
  	Timers& timers)
 :
 	extendedDomainBox_(extDomain),
+	updateInterval_(dict.getValMax<uint32>("updateInterval", 1u)),
 	particles_(prtcl),
 	geometry_(geom),
 	bTimer_("Boundary particles contact search", &timers),
-	ppTimer_("Internal particles contact search", &timers),
+	ppTimer_("Internal particles contact search", &timers)
-	dict_(dict)
 {
 
 }
@@ -69,6 +69,7 @@ bool pFlow::contactSearch::broadSearch
 		}
 		ppTimer_.end();
 		performedSearch_ = true;
+		lastUpdated_ = ti.currentIter();
 	}
 	else
 	{
@@ -92,6 +93,7 @@ bool pFlow::contactSearch::boundaryBroadSearch
 		bTimer_.start();
 		for(uint32 i=0u; i<6u; i++)
 		{
+			//output<<" boundarySearch "<< i <<" for iter "<< ti.iter()<<endl;
 			if(!BoundaryBroadSearch(
 				i,
 				ti,

src/Interaction/contactSearch/contactSearch/contactSearch.hpp

@@ -67,8 +67,6 @@ private:
 
 	Timer 		ppTimer_;
 
-	dictionary 	dict_;
-
 	virtual
 	bool BroadSearch(
 		const timeInfo& ti,
@@ -150,12 +148,6 @@ public:
 		return updateInterval_;
 	}
 
-	inline
-	const dictionary& dict()const
-	{
-		return dict_;
-	}
-
 	inline
 	const box& extendedDomainBox()const
 	{

src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp

@@ -44,9 +44,9 @@ pFlow::NBS::NBS
         position,
         flags,
         diam),
-    sizeRatio_(max(dict.getVal<real>("sizeRatio"), 1.0)),
+    sizeRatio_(max(dict.getVal<real>("sizeRatio"), one)),
-    cellExtent_(max(dict.getVal<real>("cellExtent"), 0.5)),
+    cellExtent_(max(dict.getVal<real>("cellExtent"), half)),
-    adjustableBox_(dict.getVal<Logical>("adjustableBox")),
+    adjustableBox_(false),//adjustableBox_(dict.getVal<Logical>("adjustableBox")),
     NBSLevel0_
 	(
         this->domainBox_,

src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp

@@ -31,7 +31,7 @@ pFlow::cellsWallLevel0::cellsWallLevel0
 	const ViewType1D<realx3, memory_space>& normals
 )
 :
-    cellExtent_( max(cellExtent, 0.5 ) ),
+    cellExtent_( max(cellExtent, half ) ),
     numElements_(numElements),
     numPoints_(numPoints),
     vertices_(vertices),

src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBS.cpp

@@ -23,7 +23,7 @@ Licence:
 #include "streams.hpp"
 
 pFlow::uint32 pFlow::mapperNBS::checkInterval_ = 1000;
-pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.1;
+pFlow::real pFlow::mapperNBS::enlargementFactor_ = 1.5;
 
 bool pFlow::mapperNBS::setSearchBox
 (

src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp

@@ -22,7 +22,6 @@ Licence:
 
 #include "mapperNBSKernels.hpp"
 
-
 void pFlow::mapperNBSKernels::findPointExtends
 (
     const deviceViewType1D<realx3>& points,
@@ -154,6 +153,7 @@ bool pFlow::mapperNBSKernels::buildLists
 )
 {
     auto aRange = flags.activeRange();
+    auto pp = points;
     if(flags.isAllActive() )
     { 
         Kokkos::parallel_for
@@ -162,7 +162,7 @@ bool pFlow::mapperNBSKernels::buildLists
             deviceRPolicyStatic(aRange.start(), aRange.end()),
             LAMBDA_HD(uint32 i)
             {               
-                auto ind = searchCell.pointIndex(points[i]);	
+                auto ind = searchCell.pointIndex(pp[i]);	
                 uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
                 next[i] = old;
             }

src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.hpp

@@ -152,22 +152,24 @@ public:
 		return false;
 	}
 
-
-
 	bool hearChanges
 	(
-		real t,
+		const timeInfo& ti,
-		real dt,
-		uint32 iter,
 		const message& msg, 
     	const anyList& varList
 	) override
     {
+		if(msg.equivalentTo(message::BNDR_RESET))
+		{
+			return true;
+		}
+		fatalErrorInFunction<<"Event "<< msg.eventNames() << " is not handled !"<<endl;
+		return false;
+	}
 
-		pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<<
+	bool isActive()const override
-		 msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;;
+	{
-		//notImplementedFunction;
+		return false;
-		return true;
 	}
 
 	static

src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.cpp

@@ -7,11 +7,11 @@ pFlow::boundaryGrainInteractionList<CFModel, gMModel>::boundaryGrainInteractionL
     const gMModel &geomMotion
 )
 :
-	ListPtr<boundaryGrainInteraction<CFModel,gMModel>>(6),
+	boundaryListPtr<boundaryGrainInteraction<CFModel,gMModel>>(),
 	boundaries_(grnPrtcls.pStruct().boundaries())
 {
 	//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
-	for(uint32 i=0; i<6; i++)
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,

src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.hpp

@@ -3,7 +3,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryGrainInteraction.hpp"
 
 
@@ -14,7 +14,7 @@ namespace pFlow
 template<typename contactForceModel,typename geometryMotionModel>
 class boundaryGrainInteractionList
 :
-    public ListPtr<boundaryGrainInteraction<contactForceModel,geometryMotionModel>>
+    public boundaryListPtr<boundaryGrainInteraction<contactForceModel,geometryMotionModel>>
 {
 private:
 

src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.hpp

@@ -78,14 +78,16 @@ public:
     
     ~periodicBoundaryGrainInteraction()override = default;
 
-
-    
-
     bool grainGrainInteraction(
         real dt,
 		const ContactForceModel& cfModel,
         uint32 step)override;
     
+    bool isActive()const override
+	{
+		return true;
+	}
+	
 };
 
 }

src/Interaction/grainInteraction/grainInteraction/grainInteraction.cpp

@@ -344,16 +344,12 @@ bool pFlow::grainInteraction<cFM,gMM, cLT>::afterIteration()
 template<typename cFM,typename gMM,template <class, class, class> class cLT>
 bool pFlow::grainInteraction<cFM,gMM, cLT>::hearChanges
 (
-	real t,
+	const timeInfo& ti,
-	real dt,
-	uint32 iter,
 	const message& msg, 
 	const anyList& varList
 )
 {
-	if(msg.equivalentTo(message::ITEM_REARRANGE))
+	fatalErrorInFunction<<"Event "<< msg.eventNames()<<
-	{
+	" is not handled in grainInteraction"<<endl;
-		notImplementedFunction;
+	return false;
-	}
-	return true;
 }

src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp

@@ -152,9 +152,7 @@ public:
 
 	/// Check for changes in the point structures. (overriden from observer)
 	bool hearChanges(
-		real t,
+		const timeInfo& ti,
-		real dt,
-		uint32 iter,
 		const message& msg, 
 		const anyList& varList)override;
 	

src/Interaction/grainInteraction/grainInteractionsLinearModels.cpp

@@ -60,3 +60,9 @@ createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinea
 
 createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::vibratingMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<true>>, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<true>>, pFlow::conveyorBeltMotionGeometry);
+
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGAbsoluteLinear<false>>, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::grainRolling<pFlow::cfModels::cGRelativeLinear<false>>, pFlow::conveyorBeltMotionGeometry);

src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp

@@ -41,3 +41,9 @@ Licence:
 	createBoundaryGrainInteraction(ForceModel, GeomModel)
 
 
+createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
+
+createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
+
+createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);
+createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);

src/Interaction/interaction/interaction.hpp

@@ -51,7 +51,7 @@ private:
 	const geometry& 	geometry_;
 
 	static inline
-	const message msg_ = message::ITEM_REARRANGE;
+	const message msg_ = message::ITEMS_REARRANGE;
 
 public:
 

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

@@ -148,7 +148,7 @@ public:
 		const ContactForceModel& cfModel,
 		uint32 step)
 	{
-		// for default boundary, no thing to be done 
+		// for default boundary, nothing to be done 
 		return false;
 	}
 
@@ -156,20 +156,26 @@ public:
 
 	bool hearChanges
 	(
-		real t,
+		const timeInfo& ti,
-		real dt,
-		uint32 iter,
 		const message& msg, 
     	const anyList& varList
 	) override
     {
-		
+		if(msg.equivalentTo(message::BNDR_RESET))
+		{
+			// do nothing 
+			return true;
+		}
 		pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<<
-		 msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;;
+		 msg <<endl<<" name "<< this->boundaryName() <<" type "<< this->type()<<endl;
-		//notImplementedFunction;
 		return true;
 	}
 
+	bool isActive()const override
+	{
+		return false;
+	}
+
 	static
 	uniquePtr<BoundarySphereInteractionType> create(
 		const boundaryBase& boundary,

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.cpp

@@ -7,11 +7,11 @@ pFlow::boundarySphereInteractionList<CFModel, gMModel>::boundarySphereInteractio
     const gMModel &geomMotion
 )
 :
-	ListPtr<boundarySphereInteraction<CFModel,gMModel>>(6),
+	boundaryListPtr<boundarySphereInteraction<CFModel,gMModel>>(),
 	boundaries_(sphPrtcls.pStruct().boundaries())
 {
-	//gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank());
+	output<<boundaries_.size()<<endl;
-	for(uint32 i=0; i<6; i++)
+	ForAllBoundariesPtr(i, this)
 	{
 		this->set(
 			i,

src/Interaction/sphereInteraction/boundaries/boundarySphereInteractionList.hpp

@@ -3,7 +3,7 @@
 
 
 #include "boundaryList.hpp"
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundarySphereInteraction.hpp"
 
 
@@ -14,7 +14,7 @@ namespace pFlow
 template<typename contactForceModel,typename geometryMotionModel>
 class boundarySphereInteractionList
 :
-    public ListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
+    public boundaryListPtr<boundarySphereInteraction<contactForceModel,geometryMotionModel>>
 {
 private:
 

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

@@ -78,14 +78,16 @@ public:
     
     ~periodicBoundarySphereInteraction()override = default;  
 
-
-    
-
     bool sphereSphereInteraction(
         real dt,
 		const ContactForceModel& cfModel,
         uint32 step)override;
     
+    bool isActive()const override
+	{
+		return true;
+	}
+	
 };
 
 }

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp

@@ -193,7 +193,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 	{
 		contactListMangementBoundaryTimer_.start();
 		ComputationTimer().start();
-		for(uint32 i=0; i<6u; i++)
+
+		ForAllActiveBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
 			if(activeBoundaries_[i])
 			{
@@ -202,6 +204,7 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 				BI.pwPairs().beforeBroadSearch();
 			}
 		}
+
 		ComputationTimer().end();
 		contactListMangementBoundaryTimer_.pause();	
 	}
@@ -219,7 +222,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 		fatalExit;
 	}
 	
-	for(uint32 i=0; i<6u; i++)
+	ForAllActiveBoundaries(i, boundaryInteraction_)
+	//for(size_t i=0; i<6; i++)
 	{
 		if(activeBoundaries_[i])
 		{
@@ -253,7 +257,9 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 	{
 		contactListMangementBoundaryTimer_.resume();
 		ComputationTimer().start();
-		for(uint32 i=0; i<6u; i++)
+		
+		ForAllActiveBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
 			if(activeBoundaries_[i])
 			{
@@ -262,6 +268,7 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 				BI.pwPairs().afterBroadSearch();
 			}
 		}
+
 		ComputationTimer().end();
 		contactListMangementBoundaryTimer_.end();
 	}
@@ -274,7 +281,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 	ComputationTimer().start();
 	while(requireStep.anyElement(true) && step <= 10)
 	{
-		for(uint32 i=0; i<6u; i++)
+		ForAllBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
 			if(requireStep[i] )
 			{
@@ -313,7 +321,8 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::iterate()
 	const auto& cfModel = this->forceModel_();
 	while( requireStep.anyElement(true) && step < 20 )
 	{
-		for(uint32 i=0; i<6u; i++)
+		ForAllBoundaries(i, boundaryInteraction_)
+		//for(size_t i=0; i<6; i++)
 		{
 			if(requireStep[i])
 			{
@@ -342,16 +351,18 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::afterIteration()
 template<typename cFM,typename gMM,template <class, class, class> class cLT>
 bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
 (
-	real t,
+	const timeInfo& ti,
-	real dt,
-	uint32 iter,
 	const message& msg, 
 	const anyList& varList
 )
 {
-	if(msg.equivalentTo(message::ITEM_REARRANGE))
+	if(msg.equivalentTo(message::ITEMS_REARRANGE))
 	{
 		notImplementedFunction;
+		return false;
 	}
-	return true;
+	
+	fatalErrorInFunction<<"Event "<< msg.eventNames()<<
+	" is not handled in sphereInteraction"<<endl;
+	return false;
 }

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.hpp

@@ -152,9 +152,7 @@ public:
 
 	/// Check for changes in the point structures. (overriden from observer)
 	bool hearChanges(
-		real t,
+		const timeInfo& ti,
-		real dt,
-		uint32 iter,
 		const message& msg, 
 		const anyList& varList)override;
 	

src/Interaction/sphereInteraction/sphereInteractionsLinearModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::rotation
 createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltGeometry
+createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
 //createInteraction(pFlow::cfModels::limitedLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
 //createInteraction(pFlow::cfModels::nonLimitedLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::ro
 createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
 //createInteraction(pFlow::cfModels::limitedNonLinearModNormalRolling, pFlow::multiRotationAxisMotionGeometry);
 //createInteraction(pFlow::cfModels::nonLimitedNonLinearModNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/Interaction/sphereInteraction/sphereInteractionsNonLinearModels.cpp

@@ -53,6 +53,10 @@ createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::rotat
 createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::vibratingMotionGeometry);
 createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::vibratingMotionGeometry);
 
+// conveyorBeltMotionGeometry
+createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::conveyorBeltMotionGeometry);
+createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::conveyorBeltMotionGeometry);
+
 // multiRotationAxisMotionGeometry
 //createInteraction(pFlow::cfModels::limitedNonLinearNormalRolling, pFlow::multiRotationAxisMotionGeometry);
 //createInteraction(pFlow::cfModels::nonLimitedNonLinearNormalRolling,pFlow::multiRotationAxisMotionGeometry);

src/MotionModel/CMakeLists.txt

@@ -11,6 +11,9 @@ vibratingMotion/vibratingMotion.cpp
 stationaryWall/stationaryWall.cpp
 entities/stationary/stationary.cpp
 
+conveyorBeltMotion/conveyorBeltMotion.cpp
+entities/conveyorBelt/conveyorBelt.cpp
+
 #entities/multiRotatingAxis/multiRotatingAxis.cpp
 #multiRotatingAxisMotion/multiRotatingAxisMotion.cpp
 

src/MotionModel/conveyorBeltMotion/conveyorBeltMotion.cpp

@@ -18,48 +18,58 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+#include "conveyorBeltMotion.hpp"
 
-#ifndef __eventObserver_hpp__
-#define __eventObserver_hpp__
 
-#include "eventMessage.hpp"
+pFlow::conveyorBeltMotion::conveyorBeltMotion
-
+(
-namespace pFlow
+	const objectFile &objf, 
+	repository *owner
+)
+:
+	fileDictionary(objf, owner)
 {
 
-class eventSubscriber;
+	if(!impl_readDictionary(*this) )
-
-
-class eventObserver
-{
-protected:
-	
-	const eventSubscriber* subscriber_ = nullptr;
-	
-	// if this object is linked to subscriber 
-	bool subscribed_ = false;
-
-public:
-
-	eventObserver();
-
-	eventObserver(const eventSubscriber& subscriber, bool subscribe = true );
-
-	virtual ~eventObserver();
-
-	inline bool subscribed()const {return subscribed_;}
-
-	bool subscribe(const eventSubscriber& subscriber);
-	
-	inline void invalidateSubscriber()
 	{
-		subscribed_ = false;
+		fatalErrorInFunction;
+		fatalExit;
 	}
+}
 
-	virtual bool update(const eventMessage& msg)=0;
+pFlow::conveyorBeltMotion::conveyorBeltMotion
-};
+(
+	const objectFile &objf, 
+	const dictionary &dict, 
+	repository *owner
+)
+:
+	fileDictionary(objf, dict, owner)
+{
+	if(!impl_readDictionary(*this) )
+	{
+		fatalErrorInFunction;
+		fatalExit;
+	}
+}
 
-} // pFlow
+bool pFlow::conveyorBeltMotion::write
-
+(
-
+	iOstream &os, 
-#endif // __eventObserver_hpp__
+	const IOPattern &iop
+) const
+{
+	// a global dictionary
+	dictionary newDict(fileDictionary::dictionary::name(), true);
+	if( iop.thisProcWriteData() )
+	{
+		if( !this->impl_writeDictionary(newDict) || 
+			!newDict.write(os))
+		{
+			fatalErrorInFunction<<
+			" error in writing to dictionary "<< newDict.globalName()<<endl;
+			return false;
+		}
+	}
+    return true;
+}

src/MotionModel/conveyorBeltMotion/conveyorBeltMotion.hpp

@@ -0,0 +1,104 @@
+/*------------------------------- 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 __conveyorBeltMotion_hpp__
+#define __conveyorBeltMotion_hpp__
+
+
+#include "MotionModel.hpp"
+#include "conveyorBelt.hpp"
+#include "fileDictionary.hpp"
+
+namespace pFlow
+{
+
+
+/**
+ * conveyor belt model for walls
+ * 
+ * This class is used for simulaiton that wall components are conveyor belt.
+ * 
+\verbatim
+
+// In geometryDict file, this will defines stationary walls
+...
+motionModel conveyorBelt;
+
+// this dictionary is optional
+conveyorBeltInfo
+{
+  conveyorBelt1
+  {
+    // the definition based on class conveyorBelt
+  }
+}
+...
+\endverbatim
+ *
+ */
+class conveyorBeltMotion
+:
+    public fileDictionary,
+    public MotionModel<conveyorBeltMotion, conveyorBelt>
+{
+protected:
+
+	friend MotionModel<conveyorBeltMotion, conveyorBelt>;
+  
+	bool impl_isMoving()const
+	{
+		return false;
+	}
+
+	bool impl_move(uint32, real, real)const
+	{
+		return true;
+	}
+	
+	void impl_setTime(uint32 ,real ,real )const
+	{}
+
+public:
+
+	TypeInfo("conveyorBeltMotion");
+
+	conveyorBeltMotion(const objectFile& objf, repository* owner);
+
+	conveyorBeltMotion(
+		const objectFile& objf, 
+		const dictionary& dict, 
+		repository* owner);
+
+
+    using fileDictionary::write;
+
+	bool write(iOstream& os, const IOPattern& iop)const override;
+
+	static
+    auto noneComponent()
+    {
+      return conveyorBelt();
+    }
+};
+
+} // pFlow
+
+
+#endif // __conveyorBeltMotion_hpp__

src/MotionModel/entities/conveyorBelt/conveyorBelt.cpp

@@ -18,75 +18,53 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+#include "conveyorBelt.hpp"
+#include "dictionary.hpp"
 
-#include "eventSubscriber.hpp"
-#include "Set.hpp"
 
-pFlow::eventSubscriber::~eventSubscriber()
+FUNCTION_H
+pFlow::conveyorBelt::conveyorBelt(const dictionary& dict)
 {
-	for(auto& observer:observerList_)
+  if(!read(dict))
 	{
-		observer->invalidateSubscriber();
+		fatalErrorInFunction<<
+		"  error in reading conveyorBelt from dictionary "<< dict.globalName()<<endl;
+		fatalExit;
 	}
 }
 
-bool pFlow::eventSubscriber::subscribe
+FUNCTION_H
-(
+bool pFlow::conveyorBelt::read(const dictionary& dict)
-	eventObserver* observer
-)const
 {	
-	if(observer)
+
-	{
+	tangentVelocity_ 		= dict.getVal<realx3>("tangentVelocity");
-		observerList_.push_back(observer);
+	
-		return true;
+	return true;
-	}
+}
-	else
+
+FUNCTION_H
+bool pFlow::conveyorBelt::write(dictionary& dict) const
+{
+	if( !dict.add("tangentVelocity", tangentVelocity_) )
 	{
+		fatalErrorInFunction<<
+		"  error in writing tangentVelocity to dictionary "<< dict.globalName()<<endl;
 		return false;
 	}
-}
-
-bool pFlow::eventSubscriber::unsubscribe
-(
-	eventObserver* observer
-)const
-{
-	if(observer)
-	{
-		observerList_.remove(observer);
-	}
 	return true;
 }
 
-bool pFlow::eventSubscriber::notify
+FUNCTION_H
-(
+bool pFlow::conveyorBelt::read(iIstream& is)
-	const eventMessage &msg
-)
 {
-	for ( auto& observer:observerList_ )
-	{
-		if(observer)
-			if( !observer->update(msg) ) return false;
-	}
 	
+	notImplementedFunction;
 	return true;
 }
 
-bool pFlow::eventSubscriber::notify
+FUNCTION_H
-(
+bool pFlow::conveyorBelt::write(iOstream& os)const
-	const eventMessage& msg,
-	const List<eventObserver*>& exclutionList
-)
 {
-	Set<eventObserver*> sortedExcList(exclutionList.begin(),exclutionList.end());
+  os.writeWordEntry("tangentVelocity", tangentVelocity_);
-
+  return true;
-	for(auto& observer:observerList_)
-	{
-		if( observer && sortedExcList.count(observer) == 0 )
-		{
-			if(!observer->update(msg)) return false;
-		}
-	}
-
-	return true;
 }

src/MotionModel/entities/conveyorBelt/conveyorBelt.hpp

@@ -0,0 +1,108 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __conveyorBelt_hpp__
+#define __conveyorBelt_hpp__
+
+
+#include "types.hpp"
+#include "typeInfo.hpp"
+#include "streams.hpp"
+
+namespace pFlow
+{
+
+// forward
+class dictionary;
+
+
+/**
+ * conveyor belt model for a wall
+ * 
+ */
+class conveyorBelt
+{
+private:
+
+	realx3  tangentVelocity_{0, 0, 0};
+
+public:
+
+	TypeInfoNV("conveyorBelt");
+
+	FUNCTION_HD
+	conveyorBelt()=default;
+
+	FUNCTION_H
+	explicit conveyorBelt(const dictionary& dict);
+	
+
+	FUNCTION_HD
+	conveyorBelt(const conveyorBelt&) = default;
+
+	conveyorBelt& operator=(const conveyorBelt&) = default;
+
+	INLINE_FUNCTION_HD
+	void setTime(real t)
+	{}
+
+	INLINE_FUNCTION_HD
+	realx3 linVelocityPoint(const realx3 &)const
+	{
+		return tangentVelocity_;
+	}
+
+	INLINE_FUNCTION_HD
+	realx3 transferPoint(const realx3& p, real)const
+	{
+		return p;
+	}
+
+	// - IO operation
+	FUNCTION_H 
+	bool read(const dictionary& dict);
+
+	FUNCTION_H
+	bool write(dictionary& dict) const;
+
+	FUNCTION_H
+	bool read(iIstream& is);
+
+	FUNCTION_H
+	bool write(iOstream& os)const;
+
+};
+
+inline iOstream& operator <<(iOstream& os, const conveyorBelt& obj)
+{
+	
+	return os;
+}
+
+inline iIstream& operator >>(iIstream& is, conveyorBelt& obj)
+{
+	
+	return is;
+}
+
+}
+
+
+#endif

src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp

@@ -47,7 +47,7 @@ pFlow::rotatingAxis::rotatingAxis
 	timeInterval(),
 	line(p1, p2),
 	omega_(omega),
-	rotating_(!equal(omega,0.0))
+	rotating_(!equal(omega,static_cast<real>(0.0)))
 {
 
 }
@@ -58,7 +58,7 @@ pFlow::real pFlow::rotatingAxis::setOmega(real omega)
 {
 	auto tmp = omega_;
 	omega_ = omega;
-	rotating_ = !equal(omega, 0.0);
+	rotating_ = !equal(omega, static_cast<real>(0.0));
 	return tmp;
 }
 
@@ -72,7 +72,7 @@ bool pFlow::rotatingAxis::read
 	if(!timeInterval::read(dict))return false;
 	if(!line::read(dict)) return false;
 		
-	real omega = dict.getValOrSet("omega", 0.0);
+	real omega = dict.getValOrSet("omega", static_cast<real>(0.0));
 	
 	setOmega(omega);
 	return true;

src/Particles/CMakeLists.txt

@@ -6,6 +6,9 @@ particles/shape/shape.cpp
 particles/particles.cpp
 particles/particleIdHandler/particleIdHandler.cpp
 particles/regularParticleIdHandler/regularParticleIdHandler.cpp
+
+globalDamping/globalDamping.cpp
+
 SphereParticles/sphereShape/sphereShape.cpp
 SphereParticles/sphereParticles/sphereParticles.cpp
 SphereParticles/sphereParticles/sphereParticlesKernels.cpp

src/Particles/GrainParticles/boundaryGrainParticles.hpp

@@ -50,9 +50,7 @@ public:
     const grainParticles& Particles()const;
 
     bool hearChanges(
-        real t,
+        const timeInfo& ti,
-        real dt,
-        uint32 iter,
         const message &msg,
         const anyList &varList) override
     {
@@ -65,6 +63,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+    
     static
     uniquePtr<boundaryGrainParticles> create(
         const boundaryBase &boundary,

src/Particles/GrainParticles/boundaryGrainParticlesList.cpp

@@ -5,10 +5,10 @@ pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
     grainParticles &prtcls
 )
 :
-    ListPtr(bndrs.size()), 
+    boundaryListPtr(), 
     boundaries_(bndrs)
 {
-    for(auto i=0; i<boundaries_.size(); i++)
+    ForAllBoundariesPtr(i, this)
     {
         this->set
         (
@@ -16,4 +16,5 @@ pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList(
             boundaryGrainParticles::create(boundaries_[i], prtcls)
         );
     }
+    
 }

src/Particles/GrainParticles/boundaryGrainParticlesList.hpp

@@ -3,7 +3,7 @@
 #ifndef __boundaryGrainParticlesList_hpp__
 #define __boundaryGrainParticlesList_hpp__
 
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryList.hpp"
 #include "boundaryGrainParticles.hpp"
 
@@ -12,7 +12,7 @@ namespace pFlow
 
 class boundaryGrainParticlesList
 :
-    public ListPtr<boundaryGrainParticles>
+    public boundaryListPtr<boundaryGrainParticles>
 {
 private:
 

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

@@ -142,16 +142,11 @@ pFlow::grainParticles::getParticlesInfoFromShape(
 
 pFlow::grainParticles::grainParticles(
 	systemControl &control,
-	const property& prop
+	const grainShape& gShape
 )
 :
-	particles(control),
+	particles(control, gShape),
-	grains_
+	grains_(gShape),
-	(
-		shapeFile__,
-		&control.caseSetup(),
-		prop
-	),
 	propertyId_
 	(
 		objectFile
@@ -253,7 +248,7 @@ pFlow::grainParticles::grainParticles(
 			"rVelocity",
 			dynPointStruct(),
 			intMethod,
-			rVelocity_.field()
+			rAcceleration_.field()
 		);
 
 	if( !rVelIntegration_ )
@@ -280,7 +275,7 @@ bool pFlow::grainParticles::beforeIteration()
 	particles::beforeIteration();
 	intPredictTimer_.start();
 		auto dt = this->dt();
-		dynPointStruct().predict(dt, accelertion());
+		dynPointStruct().predict(dt);
 		rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
 	intPredictTimer_.end();
 
@@ -301,8 +296,8 @@ bool pFlow::grainParticles::beforeIteration()
 bool pFlow::grainParticles::iterate() 
 {
 
-	timeInfo ti = TimeInfo();
+	const timeInfo ti = TimeInfo();
-	realx3 g = control().g();
+	const realx3 g = control().g();
 
 	particles::iterate();
 	accelerationTimer_.start();
@@ -313,25 +308,28 @@ bool pFlow::grainParticles::iterate()
 			I().deviceViewAll(),
 			contactTorque().deviceViewAll(),
 			dynPointStruct().activePointsMaskDevice(),
-			accelertion().deviceViewAll(),
+			acceleration().deviceViewAll(),
 			rAcceleration().deviceViewAll()
 			);
-		for(auto& bndry:boundaryGrainParticles_)
+		
+		ForAllActiveBoundaries(i,boundaryGrainParticles_)
 		{
-			bndry->acceleration(ti, g);
+			boundaryGrainParticles_[i].acceleration(ti, g);
 		}
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
 	
-		if(!dynPointStruct().correct(dt(), accelertion()))
+		if(!dynPointStruct().correct(ti.dt()))
 		{
 			return false;
 		}
+		real damping = dynPointStruct().dampingFactor(ti);
 		if(!rVelIntegration_().correct(
-			dt(), 
+			ti.dt(), 
 			rVelocity_, 
-			rAcceleration_))
+			rAcceleration_,
+			damping))
 		{
 			return false;
 		}

src/Particles/GrainParticles/grainParticles/grainParticles.hpp

@@ -48,7 +48,7 @@ public:
 private:
 
 	/// reference to shapes
-	ShapeType              grains_;
+	const ShapeType&       grains_;
 
 	/// property id on device
 	uint32PointField_D     propertyId_;
@@ -100,20 +100,28 @@ private:
 	  realVector&       Is,
 	  uint32Vector&     shIndex
 	);
-	/*bool initializeParticles();
 	
-	bool insertSphereParticles(
+protected:
-	    const wordVector& names,
-	    const int32IndexContainer& indices,
-	    bool setId = true);
 
-	virtual uniquePtr<List<eventObserver*>> getFieldObjectList()const override;
+	Timer& accelerationTimer()
-	*/
+	{
+		return accelerationTimer_;
+	}
+
+	Timer& intCorrectTimer()
+	{
+		return intCorrectTimer_;
+	}
+
+	integration& rVelIntegration()
+	{
+		return rVelIntegration_();
+	}
 
 public:
 
 	/// construct from systemControl and property
-	grainParticles(systemControl& control, const property& prop);
+	grainParticles(systemControl& control, const grainShape& gShape);
 
 	~grainParticles() override = default;
 
@@ -164,9 +172,7 @@ public:
 	}
 
 	bool hearChanges(
-	  real           t,
+	  const timeInfo& ti,
-	  real           dt,
-	  uint32         iter,
 	  const message& msg,
 	  const anyList& varList
 	) override
@@ -185,7 +191,12 @@ public:
 		return grainDiameter_;
 	}
 
-		const realPointField_D& coarseGrainFactor() const 
+	const realPointField_D& grainDiameter()const
+	{
+		return grainDiameter_;
+	}
+
+	const realPointField_D& coarseGrainFactor() const 
 	{
 		return coarseGrainFactor_;
 	}

src/Particles/Insertion/insertion/insertion.cpp

@@ -81,7 +81,7 @@ pFlow::insertion::readInsertionDict()
 
 	if (active_)
 	{
-    	checkForCollision_ = getVal<Logical>("checkForCollision");
+    	//checkForCollision_ = getVal<Logical>("checkForCollision");
 
 		REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
 		          << " in the simulation." << END_REPORT;

src/Particles/Insertion/insertion/insertion.hpp

@@ -42,7 +42,7 @@ private:
 	Logical    active_ = "No";
 
 	/// Check for collision? It is not active now
-	Logical    checkForCollision_ = "No";
+	Logical    checkForCollision_ = "Yes";
 
 	/// if increase velocity in case particles are failed
 	/// to be inserted due to collision

src/Particles/Insertion/insertionRegion/insertionRegion.cpp

@@ -24,6 +24,7 @@ Licence:
 #include "particles.hpp"
 #include "twoPartEntry.hpp"
 #include "types.hpp"
+#include "shape.hpp"
 
 namespace pFlow
 {

src/Particles/Insertion/insertionRegion/insertionRegion.hpp

@@ -173,7 +173,7 @@ public:
 
 	inline bool insertionTime(uint32 iter, real t, real dt) const
 	{
-		return tControl_.timeEvent(iter, t, dt);
+		return tControl_.eventTime(iter, t, dt);
 	}
 
 	uint32        numberToBeInserted(uint32 iter, real t, real dt);

src/Particles/SphereParticles/boundarySphereParticles.hpp

@@ -50,9 +50,7 @@ public:
     const sphereParticles& Particles()const;
 
     bool hearChanges(
-        real t,
+        const timeInfo& ti,
-        real dt,
-        uint32 iter,
         const message &msg,
         const anyList &varList) override
     {
@@ -65,6 +63,11 @@ public:
         return true;
     }
 
+    bool isActive()const override
+    {
+        return false;
+    }
+
     static
     uniquePtr<boundarySphereParticles> create(
         const boundaryBase &boundary,

src/Particles/SphereParticles/boundarySphereParticlesList.cpp

@@ -5,10 +5,10 @@ pFlow::boundarySphereParticlesList::boundarySphereParticlesList(
     sphereParticles &prtcls
 )
 :
-    ListPtr(bndrs.size()), 
+    boundaryListPtr(), 
     boundaries_(bndrs)
 {
-    for(auto i=0; i<boundaries_.size(); i++)
+    ForAllBoundariesPtr(i, this)
     {
         this->set
         (

src/Particles/SphereParticles/boundarySphereParticlesList.hpp

@@ -3,7 +3,7 @@
 #ifndef __boundarySphereParticlesList_hpp__
 #define __boundarySphereParticlesList_hpp__
 
-#include "ListPtr.hpp"
+#include "boundaryListPtr.hpp"
 #include "boundaryList.hpp"
 #include "boundarySphereParticles.hpp"
 
@@ -12,7 +12,7 @@ namespace pFlow
 
 class boundarySphereParticlesList
 :
-    public ListPtr<boundarySphereParticles>
+    public boundaryListPtr<boundarySphereParticles>
 {
 private:
 

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -24,180 +24,6 @@ Licence:
 #include "sphereParticlesKernels.hpp"
 
 
-//#include "setFieldList.hpp"
-/*pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
-pFlow::sphereParticles::getFieldObjectList()const
-{
-	auto objListPtr = particles::getFieldObjectList();
-
-	objListPtr().push_back(
-		static_cast<eventObserver*>(&I_) );
-
-	return objListPtr;
-}
-
-bool pFlow::sphereParticles::diameterMassInertiaPropId
-(
-	const word& shName,
-	real& diam,
-	real& mass,
-	real& I,
-	int8& propIdx
-)
-{
-	uint32 idx;
-	if( !shapes_.nameToIndex(shName, idx) )
-	{
-		printKeys(fatalErrorInFunction<<
-		"  wrong shape name in the input: "<< shName<<endl<<
-		"  available shape names are: ", shapes_.names())<<endl;
-		return false;
-	}
-
-	diam = shapes_.diameter(idx);
-	word materialName = shapes_.material(idx);
-	uint32 pIdx;
-	if( !property_.nameToIndex(materialName, pIdx) )
-	{
-		fatalErrorInFunction << 
-		"  wrong material name "<< materialName <<" specified for shape "<< shName<<
-		"  in the sphereShape dictionary.\n"<<
-		"  available materials are "<< property_.materials()<<endl;
-		return false;
-	}
-	real rho = property_.density(pIdx);
-
-	mass = Pi/6.0*pow(diam,3.0)*rho;
-	I = 0.4 * mass * pow(diam/2.0,2.0);
-	propIdx= static_cast<int8>(pIdx);
-	return true;
-}
-
-bool pFlow::sphereParticles::initializeParticles()
-{
-	
-	int32IndexContainer indices(
-		0, 
-		static_cast<int32>(shapeName_.size()));
-	
-	return insertSphereParticles(shapeName_, indices, false);
-}*/
-
-
-/*bool pFlow::sphereParticles::beforeIteration() 
-{
-	particles::beforeIteration();
-	
-	intPredictTimer_.start();
-	
-	//INFO<<"before dyn predict"<<endINFO;
-		dynPointStruct_.predict(this->dt(), accelertion_);
-	//INFO<<"after dyn predict"<<endINFO;
-
-	//INFO<<"before revel predict"<<endINFO;
-		rVelIntegration_().predict(this->dt(),rVelocity_, rAcceleration_);
-	//INFO<<"after rvel predict"<<endINFO;
-
-	intPredictTimer_.end();
-
-	return true;
-}*/
-
-
-
-
-
-/*bool pFlow::sphereParticles::afterIteration() 
-{
-	return true;
-}*/
-
-
-/*bool pFlow::sphereParticles::insertSphereParticles(
-	const wordVector& names,
-	const int32IndexContainer& indices,
-	bool setId
-	)
-{
-
-	if(names.size()!= indices.size())
-	{
-		fatalErrorInFunction <<
-		"sizes of names ("<<names.size()<<") and indices ("
-			<< indices.size()<<") do not match \n";
-		return false;
-	}
-
-	auto len = names.size();
-		
-	realVector  diamVec(len, RESERVE());
-	realVector  massVec(len, RESERVE());
-	realVector  IVec(len, RESERVE());
-	int8Vector  pIdVec(len, RESERVE());
-	int32Vector IdVec(len, RESERVE());
-
-	real d, m, I;
-	int8 pId;
-
-	ForAll(i, names )
-	{
-		
-		if (diameterMassInertiaPropId(names[i], d, m, I, pId))
-		{
-			diamVec.push_back(d);
-			massVec.push_back(m);
-			IVec.push_back(I);
-			pIdVec.push_back(pId);
-			if(setId) IdVec.push_back(idHandler_.getNextId());
-			//output<<" we are in sphereParticles nextId "<< idHandler_.nextId()<<endl;
-		}
-		else
-		{
-			fatalErrorInFunction<< "failed to calculate properties of shape " <<
-			 names[i]<<endl;
-			return false;
-		}
-		
-	}
-
-	if(!diameter_.insertSetElement(indices, diamVec))
-	{
-		fatalErrorInFunction<< " failed to insert diameters to the diameter field. \n";
-		return false;
-	}
-
-	if(!mass_.insertSetElement(indices, massVec))
-	{
-		fatalErrorInFunction<< " failed to insert mass to the mass field. \n";
-		return false;
-	}
-
-	if(!I_.insertSetElement(indices, IVec))
-	{
-		fatalErrorInFunction<< " failed to insert I to the I field. \n";
-		return false;
-	}
-
-	if(!propertyId_.insertSetElement(indices, pIdVec))
-	{
-		fatalErrorInFunction<< " failed to insert propertyId to the propertyId field. \n";
-		return false;
-	}
-
-	if(setId)
-	{
-		if( !id_.insertSetElement(indices, IdVec))
-		{
-			fatalErrorInFunction<< " failed to insert id to the id field. \n";
-			return false;
-		}	
-	}
-	
-
-	return true;
-
-}*/
-
 bool pFlow::sphereParticles::initializeParticles()
 {
 	
@@ -307,16 +133,11 @@ pFlow::sphereParticles::getParticlesInfoFromShape(
 
 pFlow::sphereParticles::sphereParticles(
 	systemControl &control,
-	const property& prop
+	const sphereShape& shpShape
 )
 :
-	particles(control),
+	particles(control, shpShape),
-	spheres_
+	spheres_(shpShape),
-	(
-		shapeFile__,
-		&control.caseSetup(),
-		prop
-	),
 	propertyId_
 	(
 		objectFile
@@ -408,7 +229,7 @@ pFlow::sphereParticles::sphereParticles(
 			"rVelocity",
 			dynPointStruct(),
 			intMethod,
-			rVelocity_.field()
+			rAcceleration_.field()
 		);
 
 	if( !rVelIntegration_ )
@@ -418,8 +239,6 @@ pFlow::sphereParticles::sphereParticles(
 		fatalExit;
 	}
 
-	WARNING<<"setFields for rVelIntegration_"<<END_WARNING;
-	
 	if(!initializeParticles())
 	{
 		fatalErrorInFunction;
@@ -428,96 +247,12 @@ pFlow::sphereParticles::sphereParticles(
 	
 } 
 
-/*bool pFlow::sphereParticles::update(const eventMessage& msg) 
-{
-	
-	if(rVelIntegration_->needSetInitialVals())
-	{
-		
-		
-		auto indexHD = pStruct().insertedPointIndex();
-	
-		auto n = indexHD.size();
-		auto index = indexHD.indicesHost();
-
-		realx3Vector rvel(n,RESERVE());
-		const auto hrVel = rVelocity_.hostView();
-
-		for(auto i=0; i<n; i++)
-		{
-			rvel.push_back( hrVel[index(i)]);
-		}
-		
-		rVelIntegration_->setInitialVals(indexHD, rvel);
-
-	}	
-
-	return true;
-}*/
-
-/*bool pFlow::sphereParticles::insertParticles
-(
-	const realx3Vector& position,
- 	const wordVector&  shapes,
- 	const setFieldList& setField
- )
-{
-	
-	if( position.size() != shapes.size() )
-	{
-		fatalErrorInFunction<< 
-		"  size of vectors position ("<<position.size()<<
-		") and shapes ("<<shapes.size()<<") are not the same. \n";
-		return false;
-	}
-
-
-
-	auto exclusionListAllPtr = getFieldObjectList(); 
-	
-	auto newInsertedPtr = pStruct().insertPoints( position, setField, time(), exclusionListAllPtr());
-	
-
-	if(!newInsertedPtr)
-	{
-		fatalErrorInFunction<<
-		"  error in inserting points into pStruct. \n";
-		return false;
-	}
-
-	auto& newInserted = newInsertedPtr();
-	
-	if(!shapeName_.insertSetElement(newInserted, shapes))
-	{
-		fatalErrorInFunction<<
-		"  error in inserting shapes into sphereParticles system.\n";
-		return false;
-	}	
-
-	if( !insertSphereParticles(shapes, newInserted) )
-	{
-		fatalErrorInFunction<< 
-		"error in inserting shapes into the sphereParticles. \n";
-		return false;
-	}
-	
-
-	auto activeR = this->activeRange();
-
-	REPORT(1)<< "Active range is "<<yellowText("["<<activeR.first<<", "<<activeR.second<<")")<<
-	" and number of active points is "<< cyanText(this->numActive())<<
-	" and pointStructure capacity is "<<cyanText(this->capacity())<<endREPORT;
-
-	return true;
-		
-}*/
-
 bool pFlow::sphereParticles::beforeIteration()
 {
 	particles::beforeIteration();
 	intPredictTimer_.start();
 		auto dt = this->dt();
-		dynPointStruct().predict(dt, accelertion());
+		dynPointStruct().predict(dt);
 		rVelIntegration_().predict(dt,rVelocity_, rAcceleration_);
 	intPredictTimer_.end();
 
@@ -537,8 +272,8 @@ bool pFlow::sphereParticles::beforeIteration()
 bool pFlow::sphereParticles::iterate() 
 {
 
-	timeInfo ti = TimeInfo();
+	const timeInfo ti = TimeInfo();
-	realx3 g = control().g();
+	const realx3 g = control().g();
 
 	particles::iterate();
 	accelerationTimer_.start();
@@ -549,25 +284,28 @@ bool pFlow::sphereParticles::iterate()
 			I().deviceViewAll(),
 			contactTorque().deviceViewAll(),
 			dynPointStruct().activePointsMaskDevice(),
-			accelertion().deviceViewAll(),
+			acceleration().deviceViewAll(),
 			rAcceleration().deviceViewAll()
 			);
-		for(auto& bndry:boundarySphereParticles_)
+		ForAllActiveBoundaries(i,boundarySphereParticles_)
 		{
-			bndry->acceleration(ti, g);
+			boundarySphereParticles_[i].acceleration(ti, g);
 		}
+
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
 	
-		if(!dynPointStruct().correct(dt(), accelertion()))
+		if(!dynPointStruct().correct(ti.dt()))
 		{
 			return false;
 		}
+		real damping = dynPointStruct().dampingFactor(ti);
 		if(!rVelIntegration_().correct(
-			dt(), 
+			ti.dt(), 
 			rVelocity_, 
-			rAcceleration_))
+			rAcceleration_,
+			damping))
 		{
 			return false;
 		}

src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp

@@ -48,7 +48,7 @@ public:
 private:
 
 	/// reference to shapes
-	ShapeType              spheres_;
+	const ShapeType&              spheres_;
 
 	/// property id on device
 	uint32PointField_D     propertyId_;
@@ -105,11 +105,26 @@ private:
 
 	virtual uniquePtr<List<eventObserver*>> getFieldObjectList()const override;
 	*/
+protected:
+	Timer& accelerationTimer()
+	{
+		return accelerationTimer_;
+	}
+
+	Timer& intCorrectTimer()
+	{
+		return intCorrectTimer_;
+	}
+
+	integration& rVelIntegration()
+	{
+		return rVelIntegration_();
+	}
 	
 public:
 
 	/// construct from systemControl and property
-	sphereParticles(systemControl& control, const property& prop);
+	sphereParticles(systemControl& control, const sphereShape& shpShape);
 
 	~sphereParticles() override = default;
 
@@ -160,9 +175,7 @@ public:
 	}
 
 	bool hearChanges(
-	  real           t,
+	  const timeInfo& ti,
-	  real           dt,
-	  uint32         iter,
 	  const message& msg,
 	  const anyList& varList
 	) override

src/Particles/dynamicPointStructure/dynamicPointStructure.cpp

@@ -21,12 +21,14 @@ Licence:
 #include "dynamicPointStructure.hpp"
 #include "systemControl.hpp"
 
+
 pFlow::dynamicPointStructure::dynamicPointStructure
 (
-	systemControl& control
+	systemControl& control,
+	real maxBSphere
 )
 :
-	pointStructure(control),
+	pointStructure(control, maxBSphere),
 	velocity_
 	(
 		objectFile(
@@ -38,6 +40,16 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		*this,
 		zero3
 	),
+	acceleration_(
+      objectFile(
+        "acceleration",
+        "",
+        objectFile::READ_IF_PRESENT,
+        objectFile::WRITE_ALWAYS
+      ),
+      *this,
+      zero3
+    ),
 	velocityUpdateTimer_("velocity boundary update", &timers()),
 	integrationMethod_
 	(
@@ -52,7 +64,7 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		"pStructPosition",
 		*this,
 		integrationMethod_,
-		pointPosition()
+		velocity_.field()
 	);
 
 	if( !integrationPos_ )
@@ -67,7 +79,7 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		"pStructVelocity",
 		*this,
 		integrationMethod_,
-		velocity_.field()
+		acceleration_.field()
 	);
 
 	if( !integrationVel_ )
@@ -77,6 +89,11 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		fatalExit;
 	}
 
+	if(control.settingsDict().containsDictionay("globalDamping"))
+	{
+		REPORT(1)<<"Reading globalDamping dictionary ..."<<END_REPORT;
+		velDamping_ = makeUnique<globalDamping>(control);
+	}
 }
 
 
@@ -85,6 +102,7 @@ bool pFlow::dynamicPointStructure::beforeIteration()
 	if(!pointStructure::beforeIteration())return false;
 	velocityUpdateTimer_.start();
 	velocity_.updateBoundariesSlaveToMasterIfRequested();
+	acceleration_.updateBoundariesSlaveToMasterIfRequested();
 	integrationPos_->updateBoundariesSlaveToMasterIfRequested();
 	integrationVel_->updateBoundariesSlaveToMasterIfRequested();
 	velocityUpdateTimer_.end();
@@ -95,34 +113,32 @@ bool pFlow::dynamicPointStructure::iterate()
 {
 	return pointStructure::iterate();
 	
-	/*real dt = this->dt();
-
-    auto& acc = time().lookupObject<realx3PointField_D>("acceleration");
-    return correct(dt, acc);*/
 }
 
-bool pFlow::dynamicPointStructure::predict(
+bool pFlow::dynamicPointStructure::afterIteration()
-    real dt,
+{
-    realx3PointField_D &acceleration)
+		//const auto ti = TimeInfo();
+		
+		auto succs = pointStructure::afterIteration();
+
+		return succs;
+}
+
+bool pFlow::dynamicPointStructure::predict(real dt)
 {
 
 	if(!integrationPos_().predict(dt, pointPosition(), velocity_ ))return false;
-	if(!integrationVel_().predict(dt, velocity_, acceleration))return false;
+	if(!integrationVel_().predict(dt, velocity_, acceleration_))return false;
 
 	return true;
 }
 
-bool pFlow::dynamicPointStructure::correct
+bool pFlow::dynamicPointStructure::correct(real dt)
-(
-	real dt,
-	realx3PointField_D& acceleration
-)
 {
-	//auto& pos = pStruct().pointPosition();
+	const auto& ti = TimeInfo();
 	
 	if(!integrationPos_().correctPStruct(dt, *this, velocity_) )return false;
-	
+	if(!integrationVel_().correct(dt, velocity_, acceleration_, dampingFactor(ti)))return false;
-	if(!integrationVel_().correct(dt, velocity_, acceleration))return false;
 
 	return true;	
 }

src/Particles/dynamicPointStructure/dynamicPointStructure.hpp

@@ -26,11 +26,13 @@ Licence:
 #include "pointFields.hpp"
 #include "integration.hpp"
 #include "uniquePtr.hpp"
+#include "globalDamping.hpp"
 
 namespace pFlow
 {
 
 class systemControl;
+//class globalDamping;
 
 class dynamicPointStructure
 :
@@ -40,10 +42,15 @@ private:
 
 	realx3PointField_D  	velocity_;
 
+	/// acceleration on device
+	realx3PointField_D      acceleration_;
+
 	uniquePtr<integration>  integrationPos_ = nullptr;
 
 	uniquePtr<integration>  integrationVel_ = nullptr;
 
+	uniquePtr<globalDamping> velDamping_ = nullptr;
+
 	Timer 					velocityUpdateTimer_;
 
 	/// @brief integration method for velocity and position
@@ -53,7 +60,7 @@ public:
 
 	TypeInfo("dynamicPointStructure");
 
-	explicit dynamicPointStructure(systemControl& control);
+	explicit dynamicPointStructure(systemControl& control, real maxBSphere);
 
 	dynamicPointStructure(const dynamicPointStructure& ps) = delete;
 		
@@ -82,6 +89,28 @@ public:
 		return velocity_;
 	}
 
+	inline 
+	const realx3PointField_D& acceleration()const
+	{
+		return acceleration_;
+	}
+
+	inline 
+	realx3PointField_D& acceleration()
+	{
+		return acceleration_;
+	}
+
+	inline
+	real dampingFactor(const timeInfo& ti)const
+	{
+		if(velDamping_)
+		{
+			return velDamping_().dampingFactor(ti);
+		}
+		return 1.0;
+	}
+
 	/// In the time loop before iterate
 	bool beforeIteration() override;
 
@@ -89,11 +118,13 @@ public:
 	/// when the component should evolve along time.
 	bool iterate() override;
 
+	bool afterIteration()override;
+	
 	/// prediction step (if any), is called in beforeIteration	
-	bool predict(real dt, realx3PointField_D& acceleration);
+	bool predict(real dt);
 
 	/// correction step, is called in iterate 
-	bool correct(real dt, realx3PointField_D& acceleration);
+	bool correct(real dt);
 
 };