positionParticles-ordered modified to accept cylinder&sphere region

.gitignore

@@ -36,6 +36,7 @@ build/**
 include/**
 bin/**
 lib/**
+test*/**
 **/**notnow
 # all possible time folders
 **/[0-9]

CMakeLists.txt

@@ -75,7 +75,7 @@ add_subdirectory(solvers)
 
 add_subdirectory(utilities)
 
-#add_subdirectory(test)
+add_subdirectory(test_newFeatures)
 
 
 install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"

benchmarks/rotatingDrum_1/caseSetup/interaction

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

benchmarks/rotatingDrum_1/caseSetup/particleInsertion

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

benchmarks/rotatingDrum_1/caseSetup/sphereShape

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

benchmarks/rotatingDrum_1/cleanThisCase

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

benchmarks/rotatingDrum_1/runThisCase

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

benchmarks/rotatingDrum_1/settings/geometryDict

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

benchmarks/rotatingDrum_1/settings/particlesDict

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

benchmarks/rotatingDrum_1/settings/settingsDict

@@ -0,0 +1,36 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	settingsDict;
+objectType 	dictionary;;
+
+run 			rotatingDrum_1;
+
+dt 				0.00001; 	// time step for integration (s)
+
+startTime 		0; 			// start time for simulation 
+
+endTime 		10; 		// end time for simulation 
+
+saveInterval 	0.2; 		// time interval for saving the simulation
+
+timePrecision   5;			// maximum number of digits for time folder 
+
+g 				(0 -9.8 0); // gravity vector (m/s2) 
+
+/*
+	Simulation domain 
+	every particles that goes outside this domain is deleted.
+*/
+domain 
+{
+	min (-0.21 -0.21 -0.01);
+	max ( 0.21  0.21  0.81);
+}
+
+integrationMethod 		AdamsBashforth3; 	// integration method 
+
+timersReport 			Yes;  	// report timers?
+
+timersReportInterval   	0.01;	// time interval for reporting timers

src/phasicFlow/CMakeLists.txt

@@ -44,6 +44,8 @@ repository/IOobject/IOobject.C
 repository/IOobject/IOfileHeader.C
 
 structuredData/box/box.C
+structuredData/cylinder/cylinder.C
+structuredData/sphere/sphere.C
 structuredData/iBox/iBoxs.C
 structuredData/line/line.C
 structuredData/pointStructure/pointStructure.C
@@ -54,6 +56,7 @@ structuredData/trisurfaceStructure/triSurface.C
 structuredData/trisurfaceStructure/multiTriSurface.C
 structuredData/trisurfaceStructure/stlFile.C
 structuredData/peakableRegion/sphereRegion/sphereRegion.C
+structuredData/peakableRegion/cylinderRegion/cylinderRegion.C
 structuredData/peakableRegion/boxRegion/boxRegion.C
 structuredData/peakableRegion/peakableRegion/peakableRegion.C
 structuredData/peakableRegion/peakableRegions.C

src/phasicFlow/structuredData/box/box.C

@@ -114,7 +114,7 @@ bool pFlow::box::write
 }
 
 FUNCTION_H
-pFlow::iIstream& pFlow::operator << (iIstream& is, box& b)
+pFlow::iIstream& pFlow::operator >>(iIstream& is, box& b)
 {
 	if(! b.read(is))
 	{
@@ -126,7 +126,7 @@ pFlow::iIstream& pFlow::operator << (iIstream& is, box& b)
 }
 
 FUNCTION_H
-pFlow::iOstream& pFlow::operator >> (iOstream& os, const box& b)
+pFlow::iOstream& pFlow::operator <<(iOstream& os, const box& b)
 {
 	
 	if(! b.write(os))

src/phasicFlow/structuredData/box/box.H

@@ -77,13 +77,13 @@ public:
 		}
 
 		INLINE_FUNCTION_HD
-		const realx3& minPoint()const
+		realx3 minPoint()const
 		{
 			return min_;
 		}
 
 		INLINE_FUNCTION_HD
-		const realx3& maxPoint()const
+		realx3 maxPoint()const
 		{
 			return max_;
 		}
@@ -103,10 +103,10 @@ public:
 };
 
 FUNCTION_H
-iIstream& operator << (iIstream& is, box& b);
+iIstream& operator >>(iIstream& is, box& b);
 
 FUNCTION_H
-iOstream& operator >> (iOstream& os, const box& b);
+iOstream& operator << (iOstream& os, const box& b);
 
 
 }

src/phasicFlow/structuredData/cylinder/cylinder.C

@@ -0,0 +1,187 @@
+/*------------------------------- 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 "cylinder.H"
+
+FUNCTION_H
+bool pFlow::cylinder::calculateParams()
+{
+
+	auto p1p2 = p2_ - p1_;
+
+	if( p1p2.length() > smallValue )
+	{
+		axisVector2_ 		= dot(p1p2,p1p2);
+		axisVector_ = p1p2;
+		return true;
+	}else
+	{
+		return false;
+	}
+}
+
+FUNCTION_H
+pFlow::cylinder::cylinder(
+	const realx3& p1,
+	const realx3& p2,
+	const real radius)
+:
+	p1_(p1),
+	p2_(p2),
+	radius2_(radius*radius)
+{
+	if(!calculateParams())
+	{
+		fatalErrorInFunction<<
+		"error in the input parameters for cylinder"<<endl;
+		fatalExit;
+	}
+}
+
+FUNCTION_H
+pFlow::cylinder::cylinder
+(
+	const dictionary & dict
+)
+:
+	p1_
+	(
+		dict.getVal<realx3>("p1")
+	),
+	p2_
+	(
+		dict.getVal<realx3>("p2")
+	)
+{
+	auto rad = dict.getVal<real>("radius");
+	radius2_= rad*rad;
+
+	if(!calculateParams())
+	{
+		fatalErrorInFunction<<
+		"error in the input parameters for cylinder in dictionary "<< dict.globalName()<<endl;
+		fatalExit;
+	}
+
+}
+
+FUNCTION_H
+pFlow::cylinder::cylinder
+(
+	iIstream& is
+)
+{
+	if( !read(is))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"error in reading cylinder from file. \n";
+		fatalExit;
+	}
+}
+
+
+FUNCTION_H
+bool pFlow::cylinder::read(iIstream & is)
+{
+	if(!is.nextData<realx3>("p1", p1_)) return false;
+	if(!is.nextData<realx3>("p2", p2_)) return false;
+	real rad;
+	if(!is.nextData<real>("radius", rad)) return false;
+	radius2_ =rad*rad;
+	return true;
+}
+
+FUNCTION_H
+bool pFlow::cylinder::write(iOstream& os)const
+{
+	os.writeWordEntry("p1", p1_);
+	os.writeWordEntry("p2", p2_);
+	os.writeWordEntry("radius", sqrt(radius2_));
+	return os.check(FUNCTION_NAME);
+}
+
+FUNCTION_H
+bool pFlow::cylinder::read
+(
+	const dictionary& dict
+)
+{
+	p1_ = dict.getVal<realx3>("p1");
+	p2_ = dict.getVal<realx3>("p2");
+	auto rad = dict.getVal<real>("radius");
+	radius2_ = rad*rad;
+	return true;
+}
+
+FUNCTION_H
+bool pFlow::cylinder::write
+(
+	dictionary& dict
+)const
+{
+	if(!dict.add("p1", p1_))
+	{
+		fatalErrorInFunction<<
+		"  error in writing p1 to dictionary "<<dict.globalName()<<endl;
+		return false;
+	}
+
+	if(!dict.add("p2", p2_))
+	{
+		fatalErrorInFunction<<
+		"  error in writing p2 to dictionary "<<dict.globalName()<<endl;
+		return false;
+	}
+
+	if(!dict.add("radius", sqrt(radius2_)) )
+	{
+		fatalErrorInFunction<<
+		"  error in writing radius to dictionary "<<dict.globalName()<<endl;
+		return false;
+	}
+	
+	return true;
+}
+
+FUNCTION_H
+pFlow::iIstream& pFlow::operator >>(iIstream& is, cylinder& b)
+{
+	if(! b.read(is))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"error in reading cylinder. \n";
+		fatalExit;
+	}
+	return is;
+}
+
+FUNCTION_H
+pFlow::iOstream& pFlow::operator << (iOstream& os, const cylinder& b)
+{
+	
+	if(! b.write(os))
+	{
+		ioErrorInFile(os.name(), os.lineNumber())<<
+		"error in writing cylinder. \n";
+		fatalExit;
+	}
+	return os;	
+}

src/phasicFlow/structuredData/cylinder/cylinder.H

@@ -0,0 +1,158 @@
+/*------------------------------- 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 __cylinder_H__
+#define __cylinder_H__
+
+#include "types.H"
+#include "dictionary.H"
+#include "iIstream.H"
+#include "iOstream.H"
+
+namespace pFlow
+{
+
+class cylinder
+{
+protected:
+	
+	// - begin point
+	realx3 	p1_;
+	
+	// - end point 
+	realx3 	p2_;
+
+	// - radius^2
+	real  	radius2_;
+
+	realx3  axisVector_;
+
+	real 		axisVector2_;
+
+	FUNCTION_H
+	bool calculateParams();
+
+public:
+
+	// - type info
+	TypeNameNV("cylinder");
+
+	//// - Constructors 
+		FUNCTION_H 
+		cylinder(const realx3& p1, const realx3& p2, const real radius);
+
+		FUNCTION_H
+		cylinder(const dictionary& dict);
+
+		FUNCTION_H
+		cylinder(iIstream& is);
+
+		FUNCTION_HD
+		cylinder(const cylinder&) = default;
+
+		FUNCTION_HD
+		cylinder(cylinder&&) = default;
+
+		FUNCTION_HD
+		cylinder& operator=(const cylinder&) = default;
+
+		FUNCTION_HD
+		cylinder& operator=(cylinder&&) = default;
+
+		~cylinder()=default;
+
+	//// - Methods 
+
+		INLINE_FUNCTION_HD
+		bool isInside(const realx3& point)const
+		{
+			auto p1Point = point-p1_;
+			auto H = cross(p1Point , axisVector_);
+			auto H2 = dot(H,H);
+			if( H2 < radius2_*axisVector2_)
+			{
+				real t = dot(p1Point, axisVector_)/axisVector2_;
+				if(t >= 0.0 && t <= 1.0)
+					return true;
+				else
+					return false;
+			}
+			else
+			{
+				return false;
+			}	
+
+		}
+
+		INLINE_FUNCTION_HD
+		const realx3& p1()const
+		{
+			return p1_;
+		}
+
+		INLINE_FUNCTION_HD
+		const realx3& p2()const
+		{
+			return p2_;
+		}
+
+		INLINE_FUNCTION_HD
+		realx3 minPoint()const
+		{
+			return min( p1_ - realx3(radius()), p2_ - realx3(radius())); // should be improved
+		}
+
+		INLINE_FUNCTION_HD
+		realx3 maxPoint()const
+		{
+			return max( p1_ + realx3(radius()), p2_ + realx3(radius())); // should be improved
+		}
+
+		INLINE_FUNCTION_HD
+		real radius()const
+		{
+			return sqrt(radius2_);
+		}
+
+	//// - IO operation 
+		FUNCTION_H
+		bool read(iIstream & is);
+
+		FUNCTION_H
+		bool write(iOstream& os)const;
+
+		FUNCTION_H
+		bool read(const dictionary& dict);
+
+		FUNCTION_H
+		bool write(dictionary& dict)const;
+};
+
+FUNCTION_H
+iIstream& operator >>(iIstream& is, cylinder& b);
+
+FUNCTION_H
+iOstream& operator << (iOstream& os, const cylinder& b);
+
+
+} // pFlow
+
+
+#endif // __cylinder_H__

src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.C

@@ -0,0 +1,73 @@
+/*------------------------------- 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 "cylinderRegion.H"
+
+
+pFlow::cylinderRegion::cylinderRegion
+(
+	const dictionary& dict
+)
+:
+	cylinder_(dict),
+	random_()
+{
+
+}
+
+bool pFlow::cylinderRegion::isInside
+(
+	const realx3& p
+) const
+{
+	return cylinder_.isInside(p);
+}
+
+pFlow::realx3 pFlow::cylinderRegion::peek()const
+{
+	for(int32 i=0; i<100;i++)
+	{
+		auto p = 
+		random_.randomNumber(cylinder_.minPoint(), cylinder_.maxPoint());
+		if( cylinder_.isInside(p)) return p;
+	}
+
+	fatalErrorInFunction<<
+	"cannot peek a random point from cylinderRegion. \n";
+	fatalExit;
+	return 0;
+}
+
+
+bool pFlow::cylinderRegion::read
+(
+	const dictionary& dict
+)
+{
+	return cylinder_.read(dict);
+}
+
+bool pFlow::cylinderRegion::write
+(
+	dictionary& dict
+)const
+{
+	return cylinder_.write(dict);
+}

src/phasicFlow/structuredData/peakableRegion/cylinderRegion/cylinderRegion.H

@@ -0,0 +1,62 @@
+/*------------------------------- 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 __cylinderRegion_H__
+#define __cylinderRegion_H__
+
+#include "cylinder.H"
+#include "uniformRandomReal.H"
+
+namespace pFlow
+{
+
+class cylinderRegion
+{
+protected:
+	cylinder 			  cylinder_;
+
+	mutable uniformRandomReal random_;
+public:
+
+	// - type info
+	TypeNameNV("cylinderRegion");
+
+	cylinderRegion(const dictionary& dict);
+
+	~cylinderRegion() = default;
+
+	//// - methods 
+	bool isInside(const realx3& p) const;
+
+	realx3 peek()const;
+
+
+	//// IO operation
+	bool read(const dictionary& dict);
+
+	bool write(dictionary& dict)const;
+
+
+};
+
+}
+
+#endif

src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.C

@@ -26,8 +26,7 @@ pFlow::sphereRegion::sphereRegion
 	const dictionary& dict
 )
 :
-	center_(dict.getVal<realx3>("center")),
+	sphere_(dict),
-	radius_(dict.getVal<real>("radius")),
 	random_()
 {
 
@@ -38,15 +37,17 @@ bool pFlow::sphereRegion::isInside
 	const realx3& p
 ) const
 {
-	return length(center_-p) <= radius_;
+	return sphere_.isInside(p);
 }
 
 pFlow::realx3 pFlow::sphereRegion::peek()const
 {
 	for(int32 i=0; i<100; ++i)
 	{
-		auto p = random_.randomNumber(center_ - realx3(radius_), center_ + realx3(radius_));
+		auto p = random_.randomNumber(
-		if(isInside(p)) return p;
+				sphere_.center() - realx3(sphere_.radius()),
+				sphere_.center() + realx3(sphere_.radius()));
+		if(sphere_.isInside(p)) return p;
 	}
 		
 	fatalErrorInFunction<<
@@ -62,8 +63,8 @@ bool pFlow::sphereRegion::read
 	const dictionary& dict
 )
 {
-	center_ = dict.getVal<realx3>("center");
+	sphere_.read(dict);
-	radius_ = dict.getVal<real>("radius");
+	
 	return true;
 }
 
@@ -72,17 +73,10 @@ bool pFlow::sphereRegion::write
 	dictionary& dict
 )const
 {
-	if(!dict.add("center", center_))
+	if(!sphere_.write(dict))
 	{
 		fatalErrorInFunction<<
-		"  error in writing center to dictionary "<<dict.globalName()<<endl;
+		"  error in writing sphere data to dictionary "<<dict.globalName()<<endl;
-		return false;
-	}
-
-	if(!dict.add("radius", radius_))
-	{
-		fatalErrorInFunction<<
-		"  error in writing radius to dictionary "<<dict.globalName()<<endl;
 		return false;
 	}
 

src/phasicFlow/structuredData/peakableRegion/sphereRegion/sphereRegion.H

@@ -22,6 +22,7 @@ Licence:
 #define __sphereRegion_H__
 
 #include "types.H"
+#include "sphere.H"
 #include "uniformRandomReal.H"
 
 namespace pFlow
@@ -33,9 +34,7 @@ class sphereRegion
 {
 protected:
 	
-	realx3 			  center_;
+	sphere 		sphere_;
-
-	real 			  radius_;
 
 	mutable uniformRandomReal random_;
 public:

src/phasicFlow/structuredData/sphere/sphere.C

@@ -0,0 +1,143 @@
+/*------------------------------- 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 "sphere.H"
+
+
+FUNCTION_H
+pFlow::sphere::sphere(
+	const realx3& center,
+	const real radius)
+:
+	center_(center),
+	radius2_(radius*radius)
+{
+	
+}
+
+FUNCTION_H
+pFlow::sphere::sphere
+(
+	const dictionary & dict
+)
+:
+	center_
+	(
+		dict.getVal<realx3>("center")
+	)
+{
+	auto rad = dict.getVal<real>("radius");
+	radius2_= rad*rad;
+}
+
+FUNCTION_H
+pFlow::sphere::sphere
+(
+	iIstream& is
+)
+{
+	if( !read(is))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"error in reading sphere from file. \n";
+		fatalExit;
+	}
+}
+
+
+FUNCTION_H
+bool pFlow::sphere::read(iIstream & is)
+{
+	if(!is.nextData<realx3>("center", center_)) return false;
+	real rad;
+	if(!is.nextData<real>("radius", rad)) return false;
+	radius2_ =rad*rad;
+	return true;
+}
+
+FUNCTION_H
+bool pFlow::sphere::write(iOstream& os)const
+{
+	os.writeWordEntry("center", center_);
+	os.writeWordEntry("radius", sqrt(radius2_));
+	return os.check(FUNCTION_NAME);
+}
+
+FUNCTION_H
+bool pFlow::sphere::read
+(
+	const dictionary& dict
+)
+{
+	center_ = dict.getVal<realx3>("center");
+	auto rad = dict.getVal<real>("radius");
+	radius2_ = rad*rad;
+	return true;
+}
+
+FUNCTION_H
+bool pFlow::sphere::write
+(
+	dictionary& dict
+)const
+{
+	if(!dict.add("center", center_))
+	{
+		fatalErrorInFunction<<
+		"  error in writing center to dictionary "<<dict.globalName()<<endl;
+		return false;
+	}
+
+
+	if(!dict.add("radius", sqrt(radius2_)) )
+	{
+		fatalErrorInFunction<<
+		"  error in writing radius to dictionary "<<dict.globalName()<<endl;
+		return false;
+	}
+	
+	return true;
+}
+
+FUNCTION_H
+pFlow::iIstream& pFlow::operator >>(iIstream& is, sphere& b)
+{
+	if(! b.read(is))
+	{
+		ioErrorInFile(is.name(), is.lineNumber())<<
+		"error in reading sphere. \n";
+		fatalExit;
+	}
+	return is;
+}
+
+FUNCTION_H
+pFlow::iOstream& pFlow::operator << (iOstream& os, const sphere& b)
+{
+	
+	if(! b.write(os))
+	{
+		ioErrorInFile(os.name(), os.lineNumber())<<
+		"error in writing sphere. \n";
+		fatalExit;
+	}
+	return os;	
+}

src/phasicFlow/structuredData/sphere/sphere.H

@@ -0,0 +1,132 @@
+/*------------------------------- 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 __sphere_H__
+#define __sphere_H__
+
+#include "types.H"
+#include "dictionary.H"
+#include "iIstream.H"
+#include "iOstream.H"
+
+namespace pFlow
+{
+
+class sphere
+{
+protected:
+	
+	// - center
+	realx3 	center_;
+	
+	// - radius^2
+	real  	radius2_;
+
+
+public:
+
+	// - type info
+	TypeNameNV("sphere");
+
+	//// - Constructors 
+		FUNCTION_H 
+		sphere(const realx3& center, const real radius);
+
+		FUNCTION_H
+		sphere(const dictionary& dict);
+
+		FUNCTION_H
+		sphere(iIstream& is);
+
+		FUNCTION_HD
+		sphere(const sphere&) = default;
+
+		FUNCTION_HD
+		sphere(sphere&&) = default;
+
+		FUNCTION_HD
+		sphere& operator=(const sphere&) = default;
+
+		FUNCTION_HD
+		sphere& operator=(sphere&&) = default;
+
+		~sphere()=default;
+
+	//// - Methods 
+
+		INLINE_FUNCTION_HD
+		bool isInside(const realx3& point)const
+		{
+			auto cPoint = point-center_;
+			auto dist2 = dot(cPoint,cPoint);
+			return dist2 < radius2_;
+
+		}
+
+		INLINE_FUNCTION_HD
+		const realx3& center()const
+		{
+			return center_;
+		}
+
+		
+		INLINE_FUNCTION_HD
+		realx3 minPoint()const
+		{
+			return center_ - realx3(radius());
+		}
+
+		INLINE_FUNCTION_HD
+		realx3 maxPoint()const
+		{
+			return center_ + realx3(radius());
+		}
+
+		INLINE_FUNCTION_HD
+		real radius()const
+		{
+			return sqrt(radius2_);
+		}
+
+	//// - IO operation 
+		FUNCTION_H
+		bool read(iIstream & is);
+
+		FUNCTION_H
+		bool write(iOstream& os)const;
+
+		FUNCTION_H
+		bool read(const dictionary& dict);
+
+		FUNCTION_H
+		bool write(dictionary& dict)const;
+};
+
+FUNCTION_H
+iIstream& operator >>(iIstream& is, sphere& b);
+
+FUNCTION_H
+iOstream& operator << (iOstream& os, const sphere& b);
+
+
+} // pFlow
+
+
+#endif // __sphere_H__

utilities/pFlowToVTK/pointFieldToVTK.H

@@ -324,7 +324,8 @@ bool convertTimeFolderPointFields(
 	auto posVec = std::as_const(pStruct).pointPosition().hostVectorAll();
 	auto* pos = posVec.data();
 
-	Report(1)<<"Writing pointStructure to vtk file."<<endReport;
+	Report(1)<<"Writing pointStructure to vtk file with "<< yellowText(pStruct.numActive())
+					 <<" active particles"<<endReport;
 	addUndstrcuturedGridField(
 		vtk(),
 		pStruct.numActive(),

utilities/particlesPhasicFlow/positionOrdered/positionOrdered.C

@@ -82,15 +82,19 @@ bool pFlow::positionOrdered::positionPointsOrdered()
 	position_.clear();
 
 	realx3 dl(diameter_);
-	auto minP = static_cast<real>(0.5)*dl + box_.minPoint();
+	auto minP =  region_->minPoint();
-	auto maxP = static_cast<real>(0.5)*dl + box_.maxPoint();
+	auto maxP =  region_->maxPoint();
 
 	auto cntr = minP;
 
 	size_t n = 0;
 	while( n < numPoints_ )
+	{
+		if(region_->isInside(cntr))
 		{
 			position_.push_back(cntr);
+			n++;
+		}
 
 		cntr += dl*uVector1_;
 		
@@ -111,7 +115,7 @@ bool pFlow::positionOrdered::positionPointsOrdered()
 				}
 			}
 		}
-		n++;
+		
 	}
 
 	return true;
@@ -139,10 +143,6 @@ pFlow::positionOrdered::positionOrdered
 	(
 		poDict_.getValOrSet("axisOrder", wordList{"x", "y", "z"})
 	),
-	box_
-	(
-		poDict_.subDict("box")
-	),
 	position_
 	(
 		maxNumberOfParticles_, RESERVE()
@@ -154,6 +154,13 @@ pFlow::positionOrdered::positionOrdered
 		fatalExit;
 	}
 
+	if(!region_)
+	{
+		fatalErrorInFunction<<"You must provided a region (box, cylinder, ...) for positioning particles in dictionary "<<
+		dict.globalName()<<endl;
+		fatalExit;
+	}
+
 	if(!positionPointsOrdered())
 	{
 		fatalExit;

utilities/particlesPhasicFlow/positionOrdered/positionOrdered.H

@@ -22,7 +22,6 @@ Licence:
 #define __positionOrdered_H__
 
 #include "positionParticles.H"
-#include "box.H"
 
 namespace pFlow
 {
@@ -42,8 +41,6 @@ protected:
 	
 	wordList 	axisOrder_;
 
-	box 		box_;
-
 	// - unit vector of the first axis
 	realx3 		uVector1_;
 

utilities/particlesPhasicFlow/positionParticles/positionParticles.C

@@ -19,6 +19,9 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "positionParticles.H"
+#include "box.H"
+#include "cylinder.H"
+#include "sphere.H"
 #include "cells.H"
 #include "contactSearchFunctions.H"
 
@@ -79,6 +82,19 @@ pFlow::positionParticles::positionParticles
 	maxNumberOfParticles_ = dict.getValOrSet("maxNumberOfParticles", static_cast<size_t>(10000));
 	
 	mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes"));
+
+	if( dict.containsDictionay("box") )	
+	{
+		region_ = makeUnique<region<box>>(dict.subDict("box"));
+	}
+	else if(dict.containsDictionay("cylinder"))
+	{
+		region_ = makeUnique<region<cylinder>>(dict.subDict("cylinder"));
+	}
+	else if(dict.containsDictionay("sphere"))
+	{
+		region_ = makeUnique<region<sphere>>(dict.subDict("sphere"));
+	}
 }
 
 

utilities/particlesPhasicFlow/positionParticles/positionParticles.H

@@ -28,11 +28,77 @@ Licence:
 namespace pFlow
 {
 
+class regionBase
+{
+public:
+	
+	regionBase() = default;
+	
+	regionBase(const regionBase&) = default;
+
+	regionBase& operator =(const regionBase&) = default;
+	
+	virtual ~regionBase() = default;
+
+	virtual bool isInside(const realx3 point)const = 0;
+
+	virtual realx3 minPoint()const =0;
+
+	virtual realx3 maxPoint()const =0;
+
+
+};
+
+template<typename T>
+class region
+:
+	public regionBase
+{
+	protected:
+
+		T 		region_;
+
+	public:
+		
+		region(const T& rgn)
+		:
+			region_(rgn)
+		{}
+
+		region(const dictionary& dict)
+		:
+			region_(dict)
+		{}
+
+		region(const region&) = default;
+		
+		region& operator =(const region&) = default;
+
+		virtual ~region()=default;
+
+		bool isInside(const realx3 point) const override
+		{
+			return region_.isInside(point);
+		}
+
+		realx3 minPoint()const override
+		{
+			return region_.minPoint();
+		}
+
+		realx3 maxPoint()const override
+		{
+			return region_.maxPoint();
+		}
+
+};
 
 class positionParticles
 {
 protected:
 
+	uniquePtr<regionBase> 		region_ = nullptr;
+
 	size_t  		maxNumberOfParticles_ = 10000;
 
 	Logical 		mortonSorting_;