correction in the macros to be compatible with OpenFOAM

DEMSystems/CMakeLists.txt

@@ -1,6 +1,8 @@
 
 set(SourceFiles 
+DEMSystem/DEMSystem.cpp
 sphereDEMSystem/sphereDEMSystem.cpp
+domainDistribute/domainDistribute.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)

DEMSystems/DEMSystem/DEMSystem.cpp

@@ -0,0 +1,75 @@
+/*------------------------------- 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 "DEMSystem.hpp"
+
+
+pFlow::coupling::DEMSystem::DEMSystem(
+	word  demSystemName,
+	int32 numDomains, 
+	const std::vector<box>& domains, 
+	int argc, 
+	char* argv[])
+:
+	ControlDict_(),
+	domains_(domains)
+{
+  
+	REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
+	Control_ = makeUnique<systemControl>(
+		ControlDict_.startTime(),
+		ControlDict_.endTime(),
+		ControlDict_.saveInterval(),
+		ControlDict_.startTimeName());
+
+}
+
+pFlow::coupling::DEMSystem::~DEMSystem()
+{}
+
+
+pFlow::uniquePtr<pFlow::coupling::DEMSystem>
+	pFlow::coupling::DEMSystem::create(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains, 
+		int argc, 
+		char* argv[]
+		)
+{
+	if( wordvCtorSelector_.search(demSystemName) )
+	{
+		return wordvCtorSelector_[demSystemName] (demSystemName, numDomains, domains, argc, argv);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< demSystemName << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		return nullptr;
+	}
+
+	return nullptr;
+}
+

DEMSystems/DEMSystem/DEMSystem.hpp

@@ -0,0 +1,128 @@
+/*------------------------------- 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 __DEMSystem_hpp__
+#define __DEMSystem_hpp__
+
+#include <vector>
+
+#include "types.hpp"
+#include "virtualConstructor.hpp"
+#include "uniquePtr.hpp"
+#include "systemControl.hpp"
+#include "readControlDict.hpp"
+
+
+namespace pFlow::coupling
+{
+
+
+class DEMSystem
+{
+protected:
+	
+	readControlDict 				ControlDict_;
+
+	uniquePtr<systemControl>		Control_ = nullptr;
+
+	std::vector<box> 				domains_;
+
+
+	// methods 
+	auto& Control()
+	{
+		return Control_();
+	}
+
+
+public:
+
+	TypeInfo("DEMSystem");
+
+	DEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains, 
+		int argc, 
+		char* argv[]);
+
+	virtual ~DEMSystem();
+
+	DEMSystem(const DEMSystem&)=delete;
+
+	DEMSystem& operator = (const DEMSystem&)=delete;
+
+	create_vCtor(
+		DEMSystem,
+		word,
+		(
+			word  demSystemName,
+			int32 numDomains, 
+			const std::vector<box>& domains, 
+			int argc, 
+			char* argv[]
+		),
+		(
+			demSystemName,
+			numDomains, 
+			domains, 
+			argc, 
+			argv
+		));
+
+	realx3 g()const
+	{
+		return Control_->g();
+	}
+
+	auto inline constexpr usingDoulle()const
+	{
+		return pFlow::usingDouble__;
+	}
+
+	virtual 
+	int32 numParInDomain(int32 di)const = 0;
+	
+	virtual
+	std::vector<int32> numParInDomain()const = 0;
+
+	virtual 
+	bool iterate(int32 n, real timeToWrite, word timeName) = 0;
+	
+	virtual 
+	real maxBounndingSphereSize()const = 0;
+
+	static 
+	uniquePtr<DEMSystem>
+		create(
+			word  demSystemName,
+			int32 numDomains, 
+			const std::vector<box>& domains, 
+			int argc, 
+			char* argv[]);
+
+
+
+};
+
+
+} // pFlow
+
+#endif // __DEMSystem_hpp__

DEMSystems/domainDistribute/domainDistribute.cpp

@@ -0,0 +1,101 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "domainDistribute.hpp"
+
+
+
+pFlow::coupling::domainDistribute::domainDistribute(
+	int32 numDomains, 
+	const Vector<box>& domains_,
+	real maxBoundingBox)
+:
+numDomains_(numDomains),
+extDomains_("extDomains", numDomains),
+particlesInDomains_("particlesInDomains", numDomains),
+numParInDomain_("numParInDomain", numDomains, 0),
+maxBoundingBoxSize_(maxBoundingBox)
+{
+
+	realx3 dl = domainExtension_ * maxBoundingBoxSize_;
+
+	for(int32 i=0; i<numDomains_; i++)
+	{
+		extDomains_[i] = extendBox(domains_[i], dl);
+	}
+}
+
+bool pFlow::coupling::domainDistribute::locateParticles(
+		ViewType1D<realx3,HostSpace> points, includeMask mask)
+	{
+		range active = mask.activeRange();
+		auto numInDomain = numParInDomain_.deviceVectorAll();
+		auto numDomains = numDomains_;
+
+		using policy = Kokkos::RangePolicy<
+			DefaultHostExecutionSpace,
+			Kokkos::IndexType<int32> >;
+
+		Kokkos::parallel_for("locateParticles",
+			policy(active.first, active.second),
+			LAMBDA_HD(int32 i){
+				if(mask(i))
+				{
+					for(int32 di=0; di<numDomains; di++)
+					{
+						if(extDomains_[i].isInside(points[i]))
+							Kokkos::atomic_add(&numInDomain[di],1);
+					}
+				}
+			});
+		Kokkos::fence();
+
+		
+		for(int32 i=0; i<numDomains_; i++)
+		{
+			particlesInDomains_[i].resize(numParInDomain_[i]);
+		}
+
+		numParInDomain_.fill(0);
+
+		auto particlesInDomainsPtr = particlesInDomains_.data();
+
+		Kokkos::parallel_for("locateParticles",
+			policy(active.first, active.second),
+			LAMBDA_HD(int32 i){
+				if(mask(i))
+				{
+					for(int32 di=0; di<numDomains; di++)
+					{
+						if(extDomains_[i].isInside(points[i]))
+						{
+
+							particlesInDomainsPtr[di][numInDomain[di]] = i;
+							Kokkos::atomic_add(&numInDomain[di],1);
+
+						}
+					}
+				}
+			});
+		Kokkos::fence();		
+
+
+		return true;
+	}

DEMSystems/domainDistribute/domainDistribute.hpp

@@ -0,0 +1,89 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __domainDistribute_hpp__
+#define __domainDistribute_hpp__
+
+#include "box.hpp"
+#include "Vectors.hpp"
+#include "VectorSingles.hpp"
+#include "pointStructure.hpp"
+
+namespace pFlow::coupling
+{
+
+class domainDistribute
+{
+protected:
+	
+	// protected members 
+
+	int32 									numDomains_;
+
+	VectorSingle<box,HostSpace> 			extDomains_;
+
+	Vector<VectorSingle<int32,HostSpace>>	particlesInDomains_;
+	
+	int32Vector_H 							numParInDomain_;
+
+
+	real maxBoundingBoxSize_;
+
+	int32 lastUpdateIter_ = 0;
+
+	int32 updateInterval_ = 1;
+
+	real domainExtension_ = 1.0;
+
+	using includeMask = typename pointStructure::activePointsHost;
+
+public:
+
+	domainDistribute(
+		int32 numDomains, 
+		const Vector<box>& domains_,
+		real maxBoundinBox);
+
+	~domainDistribute()=default;
+
+	domainDistribute(const domainDistribute&)=delete;
+
+	domainDistribute& operator=(const domainDistribute&)=delete;
+
+	int32 numDomains()const
+	{
+		return numDomains_;
+	}
+	
+	int32 numParInDomain(int32 di)const
+	{
+		return numParInDomain_[di];
+	}
+
+	//template<typename includeMask>
+	bool locateParticles(
+		ViewType1D<realx3,HostSpace> points, includeMask mask);
+	
+};
+
+
+} // pFlow::coupling
+
+#endif //__domainDistribute_hpp__

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -21,37 +21,37 @@ Licence:
 #include "sphereDEMSystem.hpp"
 
 
-pFlow::sphereDEMSystem::sphereDEMSystem(int argc, char* argv[])
+pFlow::coupling::sphereDEMSystem::sphereDEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[])
 :
-	ControlDict_()
+	DEMSystem(demSystemName, numDomains, domains, argc, argv)
 {
   
-	Report(0)<<"Initializing host/device execution spaces . . . \n";
-	Report(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endReport;
-	Report(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endReport;
+	REPORT(0)<<"Initializing host/device execution spaces . . . \n";
+	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
+	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
 
  	// initialize Kokkos
 	Kokkos::initialize( argc, argv );  
   
-	Report(0)<<"\nCreating Control repository . . ."<<endReport;
-	Control_ = makeUnique<systemControl>(
-		ControlDict_.startTime(),
-		ControlDict_.endTime(),
-		ControlDict_.saveInterval(),
-		ControlDict_.startTimeName());
+	
 
-	Report(0)<<"\nReading proprties . . . "<<endReport;
+	REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
 	property_ = makeUnique<property>(
     	Control().caseSetup().path()+propertyFile__);
 
-	Report(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endReport;
+	REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endREPORT;
 	geometry_ = geometry::create(Control(), Property());
 
-	Report(0)<<"\nReading sphere particles . . ."<<endReport;
+	REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
 	particles_ = makeUnique<sphereParticles>(Control(), Property());
 
 
-	//Report(0)<<"\nCreating particle insertion for spheres. . ."<<endReport;
+	//REPORT(0)<<"\nCreating particle insertion for spheres. . ."<<endREPORT;
 	/*insertion_ = 
 	Control().caseSetup().emplaceObject<sphereInsertion>(
 		objectFile(
@@ -64,15 +64,20 @@ pFlow::sphereDEMSystem::sphereDEMSystem(int argc, char* argv[])
 		sphParticles.shapes()
 	);*/
 
-	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 . . ."<<endREPORT;
 	interaction_ = interaction::create(
 		Control(),
 		Particles(),
 		Geometry());
 
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	particleDistribution_ = makeUnique<domainDistribute>(numDomains, domains, maxD);
+
 }
 
-pFlow::sphereDEMSystem::~sphereDEMSystem()
+pFlow::coupling::sphereDEMSystem::~sphereDEMSystem()
 {
 	interaction_.reset();
 	insertion_.reset();
@@ -85,3 +90,42 @@ pFlow::sphereDEMSystem::~sphereDEMSystem()
 	Kokkos::finalize();
 }
 
+
+
+pFlow::int32 
+	pFlow::coupling::sphereDEMSystem::numParInDomain(int32 di)const
+{
+	return particleDistribution_().numParInDomain(di);
+}
+
+std::vector<pFlow::int32> 
+	pFlow::coupling::sphereDEMSystem::numParInDomain()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;
+}
+
+
+bool pFlow::coupling::sphereDEMSystem::iterate(
+	int32 n, 
+	real timeToWrite, 
+	word timeName) 
+{
+	return true;
+}
+
+pFlow::real 
+	pFlow::coupling::sphereDEMSystem::maxBounndingSphereSize()const
+{
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	return maxD;
+}

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -21,54 +21,40 @@ Licence:
 #ifndef __sphereDEMSystem_hpp__
 #define __sphereDEMSystem_hpp__
 
-#include <array>
-
-#include "systemControl.hpp"
+#include "DEMSystem.hpp"
 #include "property.hpp"
 #include "uniquePtr.hpp"
 #include "geometry.hpp"
 #include "sphereParticles.hpp"
 #include "interaction.hpp"
 #include "Insertions.hpp"
-#include "readControlDict.hpp"
+#include "domainDistribute.hpp"
 
 
-
-
-namespace pFlow
+namespace pFlow::coupling
 {
 
 class sphereDEMSystem
+:
+	public DEMSystem
 {
 protected:
 	
-	readControlDict 						ControlDict_;
+	// protected members 
 
-	uniquePtr<systemControl>		Control_ = nullptr;
+	uniquePtr<property> 		property_ = nullptr;
 
-	uniquePtr<property> 				property_ = nullptr;
-
-	uniquePtr<geometry> 			  geometry_ = nullptr;
+	uniquePtr<geometry> 		geometry_ = nullptr;
 
 	uniquePtr<sphereParticles> 	particles_ = nullptr;
 
 	uniquePtr<sphereInsertion> 	insertion_ = nullptr;
 
-	uniquePtr<interaction> 			interaction_ = nullptr;
+	uniquePtr<interaction> 		interaction_ = nullptr;
 
+	uniquePtr<domainDistribute> particleDistribution_=nullptr;
 
-	int32 									numDomains_;
-
-	VectorDual<box> 				domains_;
-
-	Vector<int32Vector_H> 	parIndexInDomain_;
-
-	
-	auto& Control()
-	{
-		return Control_();
-	}
-
+//  protected member functions
 	auto& Property()
 	{
 		return property_();
@@ -91,27 +77,39 @@ protected:
 
 public:
 
-	sphereDEMSystem(int argc, char* argv[]);
+	TypeInfo("sphereDEMSystem");
 
-	~sphereDEMSystem();
+	sphereDEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[]);
+
+	virtual ~sphereDEMSystem();
 
 	sphereDEMSystem(const sphereDEMSystem&)=delete;
 
 	sphereDEMSystem& operator = (const sphereDEMSystem&)=delete;
 
-	
-	std::array<double,3> g()const
-	{
-		return {
-			Control_->g().x(), 
-			Control_->g().y(), 
-			Control_->g().z()};
-	}
+	add_vCtor(
+		DEMSystem,
+		sphereDEMSystem,
+		word);
 
-	bool inline usingDoulle()const
-	{
-		return pFlow::usingDouble__;
-	}
+
+	 
+	int32 numParInDomain(int32 di)const override;
+
+	
+	std::vector<int32> numParInDomain()const override;
+	
+	
+	virtual 
+	bool iterate(int32 n, real timeToWrite, word timeName) override;
+	
+	virtual 
+	real maxBounndingSphereSize()const override;
 
 };