postprocess for segregation

Update helicalMixer benchmarks

.github/scripts/sync-wiki.py

@@ -9,7 +9,7 @@ import sys
 REPO_URL = "https://github.com/PhasicFlow/phasicFlow"
 REPO_PATH = os.path.join(os.environ.get("GITHUB_WORKSPACE", ""), "repo")
 WIKI_PATH = os.path.join(os.environ.get("GITHUB_WORKSPACE", ""), "wiki")
-MAPPING_FILE = os.path.join(REPO_PATH, ".github/workflows/markdownList.yml")
+MAPPING_FILE = os.path.join(REPO_PATH, "doc/mdDocs/markdownList.yml")
 
 def load_mapping():
     """Load the markdown to wiki page mapping file."""

.github/workflows/sync-wiki.yml

@@ -7,7 +7,7 @@ on:
     paths:
       - "**/*.md"
       - ".github/workflows/sync-wiki.yml"
-      - ".github/workflows/markdownList.yml"
+      - "doc/mdDocs/markdownList.yml"
       - ".github/scripts/sync-wiki.py"
   workflow_dispatch:
 

.gitignore

@@ -37,11 +37,15 @@
 *.out
 *.app
 
-# directories
+# Exclude specific directories wherever they appear
 build/**
 include/**
 bin/**
 lib/**
+**/build/
+**/include/
+**/bin/
+**/lib/
 test*/**
 **/**notnow
 doc/code-documentation/
@@ -61,3 +65,5 @@ doc/DTAGS
 **/[0-9]*.[0-9][0-9][0-9][0-9][0-9][0-9][0-9]
 **/[0-9]*.[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]
 **/VTK
+
+

DEMSystems/DEMSystem/DEMSystem.cpp

@@ -66,12 +66,13 @@ pFlow::uniquePtr<pFlow::DEMSystem>
 		word  demSystemName,
 		const std::vector<box>& domains, 
 		int argc, 
-		char* argv[]
+		char* argv[],
+		bool requireRVel
 		)
 {
 	if( wordvCtorSelector_.search(demSystemName) )
 	{
-		return wordvCtorSelector_[demSystemName] (demSystemName, domains, argc, argv);
+		return wordvCtorSelector_[demSystemName] (demSystemName, domains, argc, argv, requireRVel);
 	}
 	else
 	{

DEMSystems/DEMSystem/DEMSystem.hpp

@@ -71,13 +71,15 @@ public:
 			word  demSystemName,
 			const std::vector<box>& domains, 
 			int argc, 
-			char* argv[]
+			char* argv[],
+			bool requireRVel
 		),
 		(
 			demSystemName,
 			domains, 
 			argc, 
-			argv
+			argv,
+			requireRVel
 		));
 
 	realx3 g()const
@@ -96,7 +98,7 @@ public:
 		return Control_();
 	}
 
-	auto inline constexpr usingDoulle()const
+	auto inline constexpr usingDouble()const
 	{
 		return pFlow::usingDouble__;
 	}
@@ -121,6 +123,9 @@ public:
 	virtual
 	span<real> diameter() = 0;
 	
+	virtual 
+	span<uint32> particleId() = 0;
+
 	virtual 
 	span<real> courseGrainFactor() = 0;
 
@@ -176,7 +181,8 @@ public:
 			word  demSystemName,
 			const std::vector<box>& domains, 
 			int argc, 
-			char* argv[]);
+			char* argv[],
+			bool requireRVel=false);
 
 };
 

DEMSystems/grainDEMSystem/grainDEMSystem.cpp

@@ -163,6 +163,12 @@ pFlow::grainDEMSystem::parIndexInDomain(int32 di)const
 	return particleDistribution_->particlesInDomain(di);
 }
 
+pFlow::span<pFlow::uint32> pFlow::grainDEMSystem::particleId()
+{
+	return span<uint32>(particleIdHost_.data(), particleIdHost_.size());
+}
+
+
 pFlow::span<pFlow::real> pFlow::grainDEMSystem::diameter() 
 {	
 	return span<real>(diameterHost_.data(), diameterHost_.size());
@@ -233,6 +239,7 @@ bool pFlow::grainDEMSystem::beforeIteration()
 	velocityHost_ = std::as_const(particles_()).velocity().hostView();
 	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
 	diameterHost_ = particles_->diameter().hostView();
+	particleIdHost_ = particles_->particleId().hostView();
 
 	if(requireRVel_)
 		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();

DEMSystems/grainDEMSystem/grainDEMSystem.hpp

@@ -63,6 +63,8 @@ protected:
 
 	ViewType1D<real, HostSpace> 	diameterHost_;
 
+	ViewType1D<uint32, HostSpace> 	particleIdHost_;
+
 	bool 							requireRVel_ = false;
 
 	ViewType1D<realx3, HostSpace> 	rVelocityHost_;
@@ -122,6 +124,8 @@ public:
 
 	span<const int32> parIndexInDomain(int32 di)const override;
 
+	span<uint32> particleId() override;
+
 	span<real> diameter() override;
 
 	span<real> courseGrainFactor() override;

DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp

@@ -165,6 +165,11 @@ pFlow::sphereDEMSystem::parIndexInDomain(int32 di)const
 	return particleDistribution_->particlesInDomain(di);
 }
 
+pFlow::span<pFlow::uint32> pFlow::sphereDEMSystem::particleId()
+{
+    return span<uint32>();
+}
+
 pFlow::span<pFlow::real> pFlow::sphereDEMSystem::diameter() 
 {	
 	return span<real>(diameterHost_.data(), diameterHost_.size());
@@ -235,6 +240,7 @@ bool pFlow::sphereDEMSystem::beforeIteration()
 	velocityHost_ = std::as_const(particles_()).velocity().hostView();
 	positionHost_ = std::as_const(particles_()).pointPosition().hostView();
 	diameterHost_ = particles_->diameter().hostView();
+	particleIdHost_ = particles_->particleId().hostView();
 
 	if(requireRVel_)
 		rVelocityHost_ = std::as_const(particles_()).rVelocity().hostView();

DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp

@@ -63,6 +63,8 @@ protected:
 
 	ViewType1D<real, HostSpace> 	diameterHost_;
 
+	ViewType1D<uint32, HostSpace> 	particleIdHost_;
+
 	bool 							requireRVel_ = false;
 
 	ViewType1D<realx3, HostSpace> 	rVelocityHost_;
@@ -122,6 +124,8 @@ public:
 
 	span<const int32> parIndexInDomain(int32 di)const override;
 
+	span<uint32> particleId() override;
+
 	span<real> diameter() override;
 
 	span<real> courseGrainFactor() override;

benchmarks/helicalMixer/caseSetup/particleInsertion

@@ -1,67 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */ 
-
-objectName particleInsertion;
-objectType dicrionary;
-
-
-active  yes;          // is insertion active?
-
-collisionCheck No;   // not implemented for yes
-
-
-
-particleInlet1
-{
-   type		      boxRegion; // type of insertion region
-   rate 	         1000000;  // insertion rate (particles/s)
-   startTime 	   0;      // (s)
-   endTime   	   2.0;    // (s)
-   interval       0.05;  //s  
-
-   boxRegionInfo 
-   {
-   	min ( -0.17 0.23 0.46); // (m,m,m)
-   	max (  0.17 0.24 0.88); // (m,m,m)
-   }
-   
-   setFields
-   {
-     velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
-   }
-   
-   mixture
-   {
-      smallParticle 1; // mixture composition of inserted particles 
-      
-   }
-}
-
-particleInlet2
-{
-   type           boxRegion; // type of insertion region
-   rate           1000000;  // insertion rate (particles/s)
-   startTime      0;      // (s)
-   endTime        2.0;    // (s)
-   interval       0.05;  //s  
-
-   boxRegionInfo 
-   {
-      min ( -0.17 0.23 0.02); // (m,m,m)
-      max (  0.17 0.24 0.44); // (m,m,m)
-   }
-   
-   setFields
-   {
-     velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
-   }
-   
-   mixture
-   {
-      largeParticle 1; // mixture composition of inserted particles 
-      
-   }
-}
-

benchmarks/helicalMixer/caseSetup/sphereShape

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

benchmarks/helicalMixer/helicalMixer_1m/caseSetup/interaction

@@ -2,11 +2,13 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-
 objectName interaction;
 objectType dicrionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
 
 materials   (glassMat wallMat); // a list of materials names
+
 densities   (2500.0 2500);      // density of materials [kg/m3]
 
 contactListType   sortedContactList;
@@ -14,6 +16,7 @@ contactListType   sortedContactList;
 model
 {
    contactForceModel    nonLinearLimited;
+   
    rollingFrictionModel normal;
 
    Yeff (1.0e6 1.0e6  // Young modulus [Pa]
@@ -28,9 +31,6 @@ model
    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);
 
@@ -42,18 +42,12 @@ model
 contactSearch
 {
     method         NBS;  
-   wallMapping    cellMapping;   
 
-   NBSInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list 
-      sizeRatio      1.1;       // bounding box size to particle diameter (max)
-   }
+    updateInterval   10;
 
-   cellMappingInfo
-   {
-      updateFrequency 10;        // each 20 timesteps, update neighbor list  
-      cellExtent     0.6;        // bounding box for particle-wall search (> 0.5)
-   }
+    sizeRatio       1.1;
 
+    cellExtent     0.55;
+
+    adjustableBox   Yes;     
 }

benchmarks/helicalMixer/helicalMixer_1m/caseSetup/particleInsertion

@@ -0,0 +1,72 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particleInsertion;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  yes;          // is insertion active?
+
+particleInlet1
+{
+   regionType     box; // type of insertion region
+   
+   rate           250000;  // insertion rate (particles/s)
+   
+   timeControl    simulationTime;
+   
+   startTime      0;      // (s)
+   
+   endTime        2.0;    // (s)
+   
+   insertionInterval    0.05;  //s  
+
+   boxInfo 
+   {
+        min (-0.17 0.23 0.46); // (m,m,m)
+        max ( 0.17 0.24 0.88); // (m,m,m)
+   }
+   
+   setFields
+   {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+   }
+   
+   mixture
+   {
+        smallParticle 1; // mixture composition of inserted particles 
+   }
+}
+
+particleInlet2
+{
+    regionType     box; // type of insertion region
+
+    rate           250000;  // insertion rate (particles/s)
+
+    timeControl    simulationTime;
+
+    startTime      0;      // (s)
+
+    endTime        2.0;    // (s)
+
+    insertionInterval    0.05;  //s  
+
+    boxInfo 
+    {
+        min ( -0.17 0.23 0.02); // (m,m,m)
+        max (  0.17 0.24 0.44); // (m,m,m)
+    }
+
+    setFields
+    {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+    }
+   
+   mixture
+   {
+        largeParticle 1; // mixture composition of inserted particles 
+   }
+}

benchmarks/helicalMixer/helicalMixer_1m/caseSetup/shapes

@@ -2,13 +2,11 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName particleInsertion;
-objectType dicrionary;
+objectName  shapes;
+objectType  dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
 
-active  no;          // is insertion active?
-
-collisionCheck No;   // not implemented for yes
-
-
+names       (smallParticle largeParticle);  // names of shapes 
+diameters   (0.004 0.00401);                // diameter of shapes 
+materials   (glassMat glassMat);            // material names for shapes 

benchmarks/helicalMixer/helicalMixer_1m/cleanThisCase

@@ -3,5 +3,5 @@ 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 
-
+rm -rf stl
 #------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_1m/runThisCase

@@ -0,0 +1,32 @@
+#!/bin/sh
+
+set -e  # Exit immediately if a command exits with a non-zero status
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "\n<--------------------------------------------------------------------->"
+echo "0) Copying stl files"
+echo "\n<--------------------------------------------------------------------->"
+mkdir -p stl
+cp -rfv $pFlow_PROJECT_DIR/resources/stls/helicalMixer/* ./stl/
+
+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
+
+echo "\n<--------------------------------------------------------------------->"
+echo "4) Converting to VtK"
+echo "<--------------------------------------------------------------------->\n"
+pFlowToVTK -f diameter id velocity --binary
+
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_1m/settings/domainDict

@@ -0,0 +1,49 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  domainDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// Simulation domain
+globalBox                                        
+{
+    min (-0.19  -0.19 -0.02);
+    max ( 0.19   0.26  0.92);
+}
+
+boundaries
+{
+    
+    left
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+}

benchmarks/helicalMixer/helicalMixer_1m/settings/geometryDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */
+objectName  geometryDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+rotatingAxisInfo
+{
+    rotAxis 
+    {
+        // end points of axis
+        p1          (0 0 0);
+        p2          (0 0 1);
+        
+        // rotation speed (rad/s) => 30 rpm
+        omega       3.1428;      
+
+        // interval for rotation of axis 
+        startTime   2.5;
+        endTime     100;
+    }
+}
+
+surfaces
+{
+    
+    helix
+    {
+        type     stlWall;       // type of the wall
+        file     helix2.stl;    // file name in stl folder
+        material wallMat;       // material name of this wall
+        motion   rotAxis;       // motion component name 
+    }
+
+    shell
+    {
+        type     stlWall;       // type of the wall
+        file     shell2.stl;    // file name in stl folder      
+        material wallMat;       // material name of this wall
+        motion   none;          // motion component name 
+    }
+        
+    plug
+    {
+        type    planeWall;
+        p1      (-0.075 -0.185 0.375);
+        p2      ( 0.075 -0.185 0.375);
+        p3      ( 0.075 -0.185 0.525);
+        p4      (-0.075 -0.185 0.525);
+        material wallMat;
+        motion   none; 
+    }
+     
+}

benchmarks/helicalMixer/helicalMixer_1m/settings/particlesDict

@@ -0,0 +1,27 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particlesDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+setFields
+{
+    defaultValue 
+    {
+        velocity        realx3  (0 0 0);        // linear velocity (m/s)
+        acceleration    realx3  (0 0 0);        // linear acceleration (m/s2)
+        rVelocity       realx3  (0 0 0);        // rotational velocity (rad/s)
+        shapeName       word    smallParticle;  // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method empty;     // creates the required fields with zero particles (empty). 
+}

benchmarks/helicalMixer/helicalMixer_1m/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  settingsDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+run         		helicalMixer;
+
+dt                  0.00001;      // time step for integration (s)
+
+startTime           0;            // start time for simulation
+
+endTime             2;          // end time for simulation
+
+saveInterval        0.05;         // time interval for saving the simulation
+
+timePrecision       4;            // maximum number of digits for time folder
+
+g                 	(0 -9.8 0);   // gravity vector (m/s2)
+
+// save necessary (i.e., required) data on disk
+includeObjects 		(diameter);                       
+
+// exclude unnecessary data from saving on disk
+excludeObjects 		(); 
+
+integrationMethod   	AdamsBashforth2; // integration method
+
+integrationHistory      off;             // Do not save integration history on the disk 
+
+writeFormat             binary;          // data writting format (ascii or binary)
+
+timersReport            Yes;             // report timers (Yes or No)
+
+timersReportInterval    0.05;            // time interval for reporting timers

benchmarks/helicalMixer/helicalMixer_250K/caseSetup/interaction

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

benchmarks/helicalMixer/helicalMixer_250K/caseSetup/particleInsertion

@@ -0,0 +1,72 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particleInsertion;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  yes;          // is insertion active?
+
+particleInlet1
+{
+   regionType     box; // type of insertion region
+   
+   rate           62500;  // insertion rate (particles/s)
+   
+   timeControl    simulationTime;
+   
+   startTime      0;      // (s)
+   
+   endTime        2.0;    // (s)
+   
+   insertionInterval    0.05;  //s  
+
+   boxInfo 
+   {
+        min (-0.17 0.23 0.46); // (m,m,m)
+        max ( 0.17 0.24 0.88); // (m,m,m)
+   }
+   
+   setFields
+   {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+   }
+   
+   mixture
+   {
+        smallParticle 1; // mixture composition of inserted particles 
+   }
+}
+
+particleInlet2
+{
+    regionType     box; // type of insertion region
+
+    rate           62500;  // insertion rate (particles/s)
+
+    timeControl    simulationTime;
+
+    startTime      0;      // (s)
+
+    endTime        2.0;    // (s)
+
+    insertionInterval    0.05;  //s  
+
+    boxInfo 
+    {
+        min ( -0.17 0.23 0.02); // (m,m,m)
+        max (  0.17 0.24 0.44); // (m,m,m)
+    }
+
+    setFields
+    {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+    }
+   
+   mixture
+   {
+        largeParticle 1; // mixture composition of inserted particles 
+   }
+}

benchmarks/helicalMixer/helicalMixer_250K/caseSetup/shapes

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

benchmarks/helicalMixer/helicalMixer_250K/cleanThisCase

@@ -0,0 +1,8 @@
+#!/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
+rm -rf stl 
+
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_250K/runThisCase

@@ -1,5 +1,12 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
+
+echo "\n<--------------------------------------------------------------------->"
+echo "0) Copying stl files"
+echo "\n<--------------------------------------------------------------------->"
+mkdir -p stl
+cp -rfv $pFlow_PROJECT_DIR/resources/stls/helicalMixer/* ./stl/
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"
 echo "<--------------------------------------------------------------------->\n"
@@ -18,6 +25,6 @@ sphereGranFlow
 echo "\n<--------------------------------------------------------------------->"
 echo "4) Converting to VtK"
 echo "<--------------------------------------------------------------------->\n"
-pFlowToVTK -f diameter id velocity
+pFlowToVTK -f diameter id velocity --binary
 
 #------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_250K/settings/domainDict

@@ -0,0 +1,49 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  domainDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// Simulation domain
+globalBox                                        
+{
+    min (-0.19  -0.19 -0.02);
+    max ( 0.19   0.26  0.92);
+}
+
+boundaries
+{
+    
+    left
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+}

benchmarks/helicalMixer/helicalMixer_250K/settings/geometryDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */
+objectName  geometryDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+rotatingAxisInfo
+{
+    rotAxis 
+    {
+        // end points of axis
+        p1          (0 0 0);
+        p2          (0 0 1);
+        
+        // rotation speed (rad/s) => 30 rpm
+        omega       3.1428;      
+
+        // interval for rotation of axis 
+        startTime   2.5;
+        endTime     100;
+    }
+}
+
+surfaces
+{
+    
+    helix
+    {
+        type     stlWall;       // type of the wall
+        file     helix2.stl;    // file name in stl folder
+        material wallMat;       // material name of this wall
+        motion   rotAxis;       // motion component name 
+    }
+
+    shell
+    {
+        type     stlWall;       // type of the wall
+        file     shell2.stl;    // file name in stl folder      
+        material wallMat;       // material name of this wall
+        motion   none;          // motion component name 
+    }
+        
+    plug
+    {
+        type    planeWall;
+        p1      (-0.075 -0.185 0.375);
+        p2      ( 0.075 -0.185 0.375);
+        p3      ( 0.075 -0.185 0.525);
+        p4      (-0.075 -0.185 0.525);
+        material wallMat;
+        motion   none; 
+    }
+     
+}

benchmarks/helicalMixer/helicalMixer_250K/settings/particlesDict

@@ -0,0 +1,27 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particlesDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+setFields
+{
+    defaultValue 
+    {
+        velocity        realx3  (0 0 0);        // linear velocity (m/s)
+        acceleration    realx3  (0 0 0);        // linear acceleration (m/s2)
+        rVelocity       realx3  (0 0 0);        // rotational velocity (rad/s)
+        shapeName       word    smallParticle;  // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method empty;     // creates the required fields with zero particles (empty). 
+}

benchmarks/helicalMixer/helicalMixer_250K/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  settingsDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+run         		helicalMixer;
+
+dt                  0.00001;      // time step for integration (s)
+
+startTime           0;            // start time for simulation
+
+endTime             7.5;          // end time for simulation
+
+saveInterval        0.05;         // time interval for saving the simulation
+
+timePrecision       4;            // maximum number of digits for time folder
+
+g                 	(0 -9.8 0);   // gravity vector (m/s2)
+
+// save necessary (i.e., required) data on disk
+includeObjects 		(diameter);                       
+
+// exclude unnecessary data from saving on disk
+excludeObjects 		(); 
+
+integrationMethod   	AdamsBashforth2; // integration method
+
+integrationHistory      off;             // Do not save integration history on the disk 
+
+writeFormat             binary;          // data writting format (ascii or binary)
+
+timersReport            Yes;             // report timers (Yes or No)
+
+timersReportInterval    0.05;            // time interval for reporting timers

benchmarks/helicalMixer/helicalMixer_2m/caseSetup/interaction

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

benchmarks/helicalMixer/helicalMixer_2m/caseSetup/particleInsertion

@@ -0,0 +1,72 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particleInsertion;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  yes;          // is insertion active?
+
+particleInlet1
+{
+   regionType     box; // type of insertion region
+   
+   rate           500000;  // insertion rate (particles/s)
+   
+   timeControl    simulationTime;
+   
+   startTime      0;      // (s)
+   
+   endTime        2.0;    // (s)
+   
+   insertionInterval    0.05;  //s  
+
+   boxInfo 
+   {
+        min (-0.17 0.23 0.46); // (m,m,m)
+        max ( 0.17 0.24 0.88); // (m,m,m)
+   }
+   
+   setFields
+   {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+   }
+   
+   mixture
+   {
+        smallParticle 1; // mixture composition of inserted particles 
+   }
+}
+
+particleInlet2
+{
+    regionType     box; // type of insertion region
+
+    rate           500000;  // insertion rate (particles/s)
+
+    timeControl    simulationTime;
+
+    startTime      0;      // (s)
+
+    endTime        2.0;    // (s)
+
+    insertionInterval    0.05;  //s  
+
+    boxInfo 
+    {
+        min ( -0.17 0.23 0.02); // (m,m,m)
+        max (  0.17 0.24 0.44); // (m,m,m)
+    }
+
+    setFields
+    {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+    }
+   
+   mixture
+   {
+        largeParticle 1; // mixture composition of inserted particles 
+   }
+}

benchmarks/helicalMixer/helicalMixer_2m/caseSetup/shapes

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

benchmarks/helicalMixer/helicalMixer_2m/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 
+rm -rf stl
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_2m/runThisCase

@@ -0,0 +1,32 @@
+#!/bin/sh
+
+set -e  # Exit immediately if a command exits with a non-zero status
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "\n<--------------------------------------------------------------------->"
+echo "0) Copying stl files"
+echo "\n<--------------------------------------------------------------------->"
+mkdir -p stl
+cp -rfv $pFlow_PROJECT_DIR/resources/stls/helicalMixer/* ./stl/
+
+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
+
+echo "\n<--------------------------------------------------------------------->"
+echo "4) Converting to VtK"
+echo "<--------------------------------------------------------------------->\n"
+pFlowToVTK -f diameter id velocity --binary
+
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_2m/settings/domainDict

@@ -0,0 +1,49 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  domainDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// Simulation domain
+globalBox                                        
+{
+    min (-0.19  -0.19 -0.02);
+    max ( 0.19   0.26  0.92);
+}
+
+boundaries
+{
+    
+    left
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+}

benchmarks/helicalMixer/helicalMixer_2m/settings/geometryDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */
+objectName  geometryDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+rotatingAxisInfo
+{
+    rotAxis 
+    {
+        // end points of axis
+        p1          (0 0 0);
+        p2          (0 0 1);
+        
+        // rotation speed (rad/s) => 30 rpm
+        omega       3.1428;      
+
+        // interval for rotation of axis 
+        startTime   2.5;
+        endTime     100;
+    }
+}
+
+surfaces
+{
+    
+    helix
+    {
+        type     stlWall;       // type of the wall
+        file     helix2.stl;    // file name in stl folder
+        material wallMat;       // material name of this wall
+        motion   rotAxis;       // motion component name 
+    }
+
+    shell
+    {
+        type     stlWall;       // type of the wall
+        file     shell2.stl;    // file name in stl folder      
+        material wallMat;       // material name of this wall
+        motion   none;          // motion component name 
+    }
+        
+    plug
+    {
+        type    planeWall;
+        p1      (-0.075 -0.185 0.375);
+        p2      ( 0.075 -0.185 0.375);
+        p3      ( 0.075 -0.185 0.525);
+        p4      (-0.075 -0.185 0.525);
+        material wallMat;
+        motion   none; 
+    }
+     
+}

benchmarks/helicalMixer/helicalMixer_2m/settings/particlesDict

@@ -0,0 +1,27 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particlesDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+setFields
+{
+    defaultValue 
+    {
+        velocity        realx3  (0 0 0);        // linear velocity (m/s)
+        acceleration    realx3  (0 0 0);        // linear acceleration (m/s2)
+        rVelocity       realx3  (0 0 0);        // rotational velocity (rad/s)
+        shapeName       word    smallParticle;  // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method empty;     // creates the required fields with zero particles (empty). 
+}

benchmarks/helicalMixer/helicalMixer_2m/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  settingsDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+run         		helicalMixer;
+
+dt                  0.00001;      // time step for integration (s)
+
+startTime           0;            // start time for simulation
+
+endTime             2;          // end time for simulation
+
+saveInterval        0.05;         // time interval for saving the simulation
+
+timePrecision       4;            // maximum number of digits for time folder
+
+g                 	(0 -9.8 0);   // gravity vector (m/s2)
+
+// save necessary (i.e., required) data on disk
+includeObjects 		(diameter);                       
+
+// exclude unnecessary data from saving on disk
+excludeObjects 		(); 
+
+integrationMethod   	AdamsBashforth2; // integration method
+
+integrationHistory      off;             // Do not save integration history on the disk 
+
+writeFormat             binary;          // data writting format (ascii or binary)
+
+timersReport            Yes;             // report timers (Yes or No)
+
+timersReportInterval    0.05;            // time interval for reporting timers

benchmarks/helicalMixer/helicalMixer_4m/caseSetup/interaction

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

benchmarks/helicalMixer/helicalMixer_4m/caseSetup/particleInsertion

@@ -0,0 +1,72 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particleInsertion;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  yes;          // is insertion active?
+
+particleInlet1
+{
+   regionType     box; // type of insertion region
+   
+   rate           1000000;  // insertion rate (particles/s)
+   
+   timeControl    simulationTime;
+   
+   startTime      0;      // (s)
+   
+   endTime        2.0;    // (s)
+   
+   insertionInterval    0.05;  //s  
+
+   boxInfo 
+   {
+        min (-0.17 0.23 0.46); // (m,m,m)
+        max ( 0.17 0.24 0.88); // (m,m,m)
+   }
+   
+   setFields
+   {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+   }
+   
+   mixture
+   {
+        smallParticle 1; // mixture composition of inserted particles 
+   }
+}
+
+particleInlet2
+{
+    regionType     box; // type of insertion region
+
+    rate           1000000;  // insertion rate (particles/s)
+
+    timeControl    simulationTime;
+
+    startTime      0;      // (s)
+
+    endTime        2.0;    // (s)
+
+    insertionInterval    0.05;  //s  
+
+    boxInfo 
+    {
+        min ( -0.17 0.23 0.02); // (m,m,m)
+        max (  0.17 0.24 0.44); // (m,m,m)
+    }
+
+    setFields
+    {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+    }
+   
+   mixture
+   {
+        largeParticle 1; // mixture composition of inserted particles 
+   }
+}

benchmarks/helicalMixer/helicalMixer_4m/caseSetup/shapes

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

benchmarks/helicalMixer/helicalMixer_4m/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 
+rm -rf stl
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_4m/runThisCase

@@ -0,0 +1,32 @@
+#!/bin/sh
+
+set -e  # Exit immediately if a command exits with a non-zero status
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "\n<--------------------------------------------------------------------->"
+echo "0) Copying stl files"
+echo "\n<--------------------------------------------------------------------->"
+mkdir -p stl
+cp -rfv $pFlow_PROJECT_DIR/resources/stls/helicalMixer/* ./stl/
+
+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
+
+echo "\n<--------------------------------------------------------------------->"
+echo "4) Converting to VtK"
+echo "<--------------------------------------------------------------------->\n"
+pFlowToVTK -f diameter id velocity --binary
+
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_4m/settings/domainDict

@@ -0,0 +1,49 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  domainDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// Simulation domain
+globalBox                                        
+{
+    min (-0.19  -0.19 -0.02);
+    max ( 0.19   0.26  0.92);
+}
+
+boundaries
+{
+    
+    left
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+}

benchmarks/helicalMixer/helicalMixer_4m/settings/geometryDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */
+objectName  geometryDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+rotatingAxisInfo
+{
+    rotAxis 
+    {
+        // end points of axis
+        p1          (0 0 0);
+        p2          (0 0 1);
+        
+        // rotation speed (rad/s) => 30 rpm
+        omega       3.1428;      
+
+        // interval for rotation of axis 
+        startTime   2.5;
+        endTime     100;
+    }
+}
+
+surfaces
+{
+    
+    helix
+    {
+        type     stlWall;       // type of the wall
+        file     helix2.stl;    // file name in stl folder
+        material wallMat;       // material name of this wall
+        motion   rotAxis;       // motion component name 
+    }
+
+    shell
+    {
+        type     stlWall;       // type of the wall
+        file     shell2.stl;    // file name in stl folder      
+        material wallMat;       // material name of this wall
+        motion   none;          // motion component name 
+    }
+        
+    plug
+    {
+        type    planeWall;
+        p1      (-0.075 -0.185 0.375);
+        p2      ( 0.075 -0.185 0.375);
+        p3      ( 0.075 -0.185 0.525);
+        p4      (-0.075 -0.185 0.525);
+        material wallMat;
+        motion   none; 
+    }
+     
+}

benchmarks/helicalMixer/helicalMixer_4m/settings/particlesDict

@@ -0,0 +1,27 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particlesDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+setFields
+{
+    defaultValue 
+    {
+        velocity        realx3  (0 0 0);        // linear velocity (m/s)
+        acceleration    realx3  (0 0 0);        // linear acceleration (m/s2)
+        rVelocity       realx3  (0 0 0);        // rotational velocity (rad/s)
+        shapeName       word    smallParticle;  // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method empty;     // creates the required fields with zero particles (empty). 
+}

benchmarks/helicalMixer/helicalMixer_4m/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  settingsDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+run         		helicalMixer;
+
+dt                  0.00001;      // time step for integration (s)
+
+startTime           0;            // start time for simulation
+
+endTime             2;          // end time for simulation
+
+saveInterval        0.05;         // time interval for saving the simulation
+
+timePrecision       4;            // maximum number of digits for time folder
+
+g                 	(0 -9.8 0);   // gravity vector (m/s2)
+
+// save necessary (i.e., required) data on disk
+includeObjects 		(diameter);                       
+
+// exclude unnecessary data from saving on disk
+excludeObjects 		(); 
+
+integrationMethod   	AdamsBashforth2; // integration method
+
+integrationHistory      off;             // Do not save integration history on the disk 
+
+writeFormat             binary;          // data writting format (ascii or binary)
+
+timersReport            Yes;             // report timers (Yes or No)
+
+timersReportInterval    0.05;            // time interval for reporting timers

benchmarks/helicalMixer/helicalMixer_500k/caseSetup/interaction

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

benchmarks/helicalMixer/helicalMixer_500k/caseSetup/particleInsertion

@@ -0,0 +1,72 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particleInsertion;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+active  yes;          // is insertion active?
+
+particleInlet1
+{
+   regionType     box; // type of insertion region
+   
+   rate           125000;  // insertion rate (particles/s)
+   
+   timeControl    simulationTime;
+   
+   startTime      0;      // (s)
+   
+   endTime        2.0;    // (s)
+   
+   insertionInterval    0.05;  //s  
+
+   boxInfo 
+   {
+        min (-0.17 0.23 0.46); // (m,m,m)
+        max ( 0.17 0.24 0.88); // (m,m,m)
+   }
+   
+   setFields
+   {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+   }
+   
+   mixture
+   {
+        smallParticle 1; // mixture composition of inserted particles 
+   }
+}
+
+particleInlet2
+{
+    regionType     box; // type of insertion region
+
+    rate           125000;  // insertion rate (particles/s)
+
+    timeControl    simulationTime;
+
+    startTime      0;      // (s)
+
+    endTime        2.0;    // (s)
+
+    insertionInterval    0.05;  //s  
+
+    boxInfo 
+    {
+        min ( -0.17 0.23 0.02); // (m,m,m)
+        max (  0.17 0.24 0.44); // (m,m,m)
+    }
+
+    setFields
+    {
+        velocity    realx3 (0.0 -0.3 0.0); // initial velocity of inserted particles 
+    }
+   
+   mixture
+   {
+        largeParticle 1; // mixture composition of inserted particles 
+   }
+}

benchmarks/helicalMixer/helicalMixer_500k/caseSetup/shapes

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

benchmarks/helicalMixer/helicalMixer_500k/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 
+rm -rf stl
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_500k/runThisCase

@@ -0,0 +1,32 @@
+#!/bin/sh
+
+set -e  # Exit immediately if a command exits with a non-zero status
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "\n<--------------------------------------------------------------------->"
+echo "0) Copying stl files"
+echo "\n<--------------------------------------------------------------------->"
+mkdir -p stl
+cp -rfv $pFlow_PROJECT_DIR/resources/stls/helicalMixer/* ./stl/
+
+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
+
+echo "\n<--------------------------------------------------------------------->"
+echo "4) Converting to VtK"
+echo "<--------------------------------------------------------------------->\n"
+pFlowToVTK -f diameter id velocity --binary
+
+#------------------------------------------------------------------------------

benchmarks/helicalMixer/helicalMixer_500k/settings/domainDict

@@ -0,0 +1,49 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  domainDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// Simulation domain
+globalBox                                        
+{
+    min (-0.19  -0.19 -0.02);
+    max ( 0.19   0.26  0.92);
+}
+
+boundaries
+{
+    
+    left
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    right
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    bottom
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    top
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    rear
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+
+    front
+    {
+        type     exit;        // other options: periodic, reflective 
+    }
+}

benchmarks/helicalMixer/helicalMixer_500k/settings/geometryDict

@@ -0,0 +1,60 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */
+objectName  geometryDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+// motion model: rotating object around an axis
+motionModel rotatingAxis; 
+
+rotatingAxisInfo
+{
+    rotAxis 
+    {
+        // end points of axis
+        p1          (0 0 0);
+        p2          (0 0 1);
+        
+        // rotation speed (rad/s) => 30 rpm
+        omega       3.1428;      
+
+        // interval for rotation of axis 
+        startTime   2.5;
+        endTime     100;
+    }
+}
+
+surfaces
+{
+    
+    helix
+    {
+        type     stlWall;       // type of the wall
+        file     helix2.stl;    // file name in stl folder
+        material wallMat;       // material name of this wall
+        motion   rotAxis;       // motion component name 
+    }
+
+    shell
+    {
+        type     stlWall;       // type of the wall
+        file     shell2.stl;    // file name in stl folder      
+        material wallMat;       // material name of this wall
+        motion   none;          // motion component name 
+    }
+        
+    plug
+    {
+        type    planeWall;
+        p1      (-0.075 -0.185 0.375);
+        p2      ( 0.075 -0.185 0.375);
+        p3      ( 0.075 -0.185 0.525);
+        p4      (-0.075 -0.185 0.525);
+        material wallMat;
+        motion   none; 
+    }
+     
+}

benchmarks/helicalMixer/helicalMixer_500k/settings/particlesDict

@@ -0,0 +1,27 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */ 
+objectName  particlesDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+setFields
+{
+    defaultValue 
+    {
+        velocity        realx3  (0 0 0);        // linear velocity (m/s)
+        acceleration    realx3  (0 0 0);        // linear acceleration (m/s2)
+        rVelocity       realx3  (0 0 0);        // rotational velocity (rad/s)
+        shapeName       word    smallParticle;  // name of the particle shape 
+    }
+
+    selectors
+    {}
+}
+
+positionParticles
+{
+    method empty;     // creates the required fields with zero particles (empty). 
+}

benchmarks/helicalMixer/helicalMixer_500k/settings/settingsDict

@@ -0,0 +1,37 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName  settingsDict;
+objectType  dictionary;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+run         		helicalMixer;
+
+dt                  0.00001;      // time step for integration (s)
+
+startTime           0;            // start time for simulation
+
+endTime             2;          // end time for simulation
+
+saveInterval        0.05;         // time interval for saving the simulation
+
+timePrecision       4;            // maximum number of digits for time folder
+
+g                 	(0 -9.8 0);   // gravity vector (m/s2)
+
+// save necessary (i.e., required) data on disk
+includeObjects 		(diameter);                       
+
+// exclude unnecessary data from saving on disk
+excludeObjects 		(); 
+
+integrationMethod   	AdamsBashforth2; // integration method
+
+integrationHistory      off;             // Do not save integration history on the disk 
+
+writeFormat             binary;          // data writting format (ascii or binary)
+
+timersReport            Yes;             // report timers (Yes or No)
+
+timersReportInterval    0.05;            // time interval for reporting timers

benchmarks/helicalMixer/settings/geometryDict

@@ -1,56 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	geometryDict;
-objectType 	dictionary;
-
-// motion model: rotating object around an axis
-motionModel rotatingAxisMotion; 
-
-surfaces
-{
-	
-	helix
-	{
-		type 	 stlWall;   // type of the wall
-		file 	 helix2.stl;	// file name in stl folder
-		material wallMat;   // material name of this wall
-		motion 	 rotAxis;	// motion component name 
-	}
-
-	
-	shell
-	{
-		type 	 stlWall;  	// type of the wall
-		file 	 shell2.stl;	// file name in stl folder		
-		material wallMat;   // material name of this wall
-		motion 	 none;		// motion component name 
-	}
-
-		
-	plug
-	{
-		type 	planeWall;
-		p1 (-0.075 -0.185 0.375);
-		p2 ( 0.075 -0.185 0.375);
-		p3 ( 0.075 -0.185 0.525);
-		p4 (-0.075 -0.185 0.525);
-		material wallMat;   // material name of this wall
-		motion 	 none;		// motion component name 
-	}
-	
-	
-}
-
-// information for rotatingAxisMotion motion model 
-rotatingAxisMotionInfo
-{
-	rotAxis 
-	{
-		p1 ( 0 0 0);
-		p2 ( 0 0 1);
-		omega 0; //3.1428; 		// rotation speed (rad/s) => 30 rpm
-	}
-}

benchmarks/helicalMixer/settings/particlesDict

@@ -1,31 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-
-objectName 	particlesDict;
-objectType 	dictionary;
-
-setFields
-{
-	defaultValue 
-	{
-		velocity 		realx3 	(0 0 0); 		// linear velocity (m/s)
-		acceleration 	realx3 	(0 0 0); 		// linear acceleration (m/s2)
-		rotVelocity 	realx3 	(0 0 0);  		// rotational velocity (rad/s)
-		shapeName 		word	smallParticle; 	// name of the particle shape 
-	}
-
-	selectors
-	{}
-}
-
-// positions particles 
-positionParticles
-{
-	method empty;     		// creates the required fields with zero particles (empty). 
-
-	maxNumberOfParticles 	4100000; // maximum number of particles in the simulation
-	mortonSorting 			Yes;    // perform initial sorting based on morton code?   
-
-}

benchmarks/helicalMixer/settings/settingsDict

@@ -1,36 +0,0 @@
-/* -------------------------------*- C++ -*--------------------------------- *\ 
-|  phasicFlow File                                                            | 
-|  copyright: www.cemf.ir                                                     | 
-\* ------------------------------------------------------------------------- */  
-objectName 	settingsDict;
-objectType 	dictionary;;
-
-run 			inclinedScrewConveyor;
-
-dt 				0.00001; 	// time step for integration (s)
-
-startTime 		2.9; 			// start time for simulation 
-
-endTime 		7;			// end time for simulation 
-
-saveInterval 	0.05; 		// time interval for saving the simulation
-
-timePrecision   3;			// 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.19  -0.19 -0.02);
-	max ( 0.19   0.26  0.92);
-}
-
-integrationMethod 		AdamsBashforth2; 	// integration method 
-
-timersReport 			Yes;  	// report timers?
-
-timersReportInterval   	0.01;	// time interval for reporting timers

benchmarks/rotatingDrum/rotatingDrum_1mParticles/runThisCase

@@ -1,4 +1,5 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"

benchmarks/rotatingDrum/rotatingDrum_250kParticles/runThisCase

@@ -1,4 +1,5 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"

benchmarks/rotatingDrum/rotatingDrum_2mParticles/runThisCase

@@ -1,4 +1,5 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"

benchmarks/rotatingDrum/rotatingDrum_4mParticles/runThisCase

@@ -1,4 +1,5 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"

benchmarks/rotatingDrum/rotatingDrum_8mParticles/runThisCase

@@ -1,4 +1,5 @@
 #!/bin/sh
+set -e  # Exit immediately if a command exits with a non-zero status
 cd ${0%/*} || exit 1    # Run from this directory
 echo "\n<--------------------------------------------------------------------->"
 echo "1) Creating particles"

cmake/bashrc

@@ -19,7 +19,7 @@ export pFlow_SRC_DIR="$pFlow_PROJECT_DIR/src"
 
 export Kokkos_DIR="$kokkosDir"
 
-export Zoltan_DIR="$projectDir/Zoltan"
+#export Zoltan_DIR="$projectDir/Zoltan"
 
 # Cleanup variables (done as final statement for a clean exit code)
 unset projectDir

cmake/zoltanInstallCheck.cmake

@@ -0,0 +1,44 @@
+# Macro to check for Zoltan installation and build it if needed
+# Usage: zoltan_find_or_build(ZOLTAN_DIR)
+# Returns: ZOLTAN_INCLUDE_DIR, ZOLTAN_LIBRARY
+
+macro(zoltan_find_or_build ZOLTAN_DIR)
+    # Set the Zoltan directory
+    set(ZOLTAN_PREFIX "${ZOLTAN_DIR}" CACHE STRING "Zoltan install directory")
+    message(STATUS "Zoltan install directory is ${ZOLTAN_PREFIX}")
+
+    # Check if the Zoltan library is already built
+    find_path(ZOLTAN_INCLUDE_DIR zoltan.h PATHS "${ZOLTAN_PREFIX}/include")
+    message(STATUS "Zoltan include path: ${ZOLTAN_INCLUDE_DIR}")
+
+    find_library(ZOLTAN_LIBRARY zoltan PATHS "${ZOLTAN_PREFIX}/lib")
+    message(STATUS "Zoltan lib path: ${ZOLTAN_LIBRARY}")
+    
+    # Check if Zoltan library exists, if not compile it using buildlib script
+    if(NOT ZOLTAN_LIBRARY)
+        message(STATUS "Zoltan library not found. Compiling from source using buildlib script...")
+
+        # Execute the buildlib bash script
+        execute_process(
+            COMMAND bash ${ZOLTAN_PREFIX}/buildlib
+            WORKING_DIRECTORY ${ZOLTAN_PREFIX}
+            RESULT_VARIABLE ZOLTAN_BUILD_RESULT
+            OUTPUT_VARIABLE ZOLTAN_BUILD_OUTPUT
+            ERROR_VARIABLE ZOLTAN_BUILD_ERROR
+        )
+
+        if(NOT ZOLTAN_BUILD_RESULT EQUAL 0)
+            message(FATAL_ERROR "Failed to build Zoltan library using buildlib script. Error: ${ZOLTAN_BUILD_ERROR}")
+        endif()
+
+        # Try to find the library again after building
+        find_library(ZOLTAN_LIBRARY zoltan PATHS "${ZOLTAN_PREFIX}/lib" NO_DEFAULT_PATH)
+        find_path(ZOLTAN_INCLUDE_DIR zoltan.h PATHS "${ZOLTAN_PREFIX}/include" NO_DEFAULT_PATH)
+
+        if(NOT ZOLTAN_LIBRARY)
+            message(FATAL_ERROR "Failed to locate Zoltan library after building")
+        endif()
+
+        message(STATUS "Successfully built Zoltan library at ${ZOLTAN_LIBRARY}")
+    endif()
+endmacro()

src/Integration/AdamsBashforth2/processorAB2BoundaryIntegration.cpp

@@ -0,0 +1,71 @@
+#include "processorAB2BoundaryIntegration.hpp"
+#include "AdamsBashforth2.hpp"
+#include "AB2Kernels.hpp"
+#include "boundaryConfigs.hpp"
+
+pFlow::processorAB2BoundaryIntegration::processorAB2BoundaryIntegration(
+    const boundaryBase &boundary, 
+    const pointStructure &pStruct,  
+    const word &method,
+    integration& intgrtn
+)
+:
+    boundaryIntegration(boundary, pStruct, method, intgrtn)
+{}
+
+bool pFlow::processorAB2BoundaryIntegration::correct(
+    real dt, 
+	const realx3PointField_D& y, 
+	const realx3PointField_D& dy
+)
+{
+
+#ifndef BoundaryModel1
+    if(this->isBoundaryMaster())
+    {
+        const uint32 thisIndex = thisBoundaryIndex();
+		const auto&  AB2 =  static_cast<const AdamsBashforth2&>(Integration());
+        const auto& dy1View = AB2.BoundaryField(thisIndex).neighborProcField().deviceView();
+		const auto& dyView = dy.BoundaryField(thisIndex).neighborProcField().deviceView();
+		const auto& yView = y.BoundaryField(thisIndex).neighborProcField().deviceView();        
+        const rangeU32 aRange(0u, dy1View.size());
+        return  AB2Kernels::intAllActive(
+			"AB2Integration::correct."+this->boundaryName(), 
+			dt,
+            aRange, 
+			yView, 
+			dyView, 
+			dy1View
+		); 
+    }
+#endif //BoundaryModel1
+    
+    
+    return true;
+}
+
+bool pFlow::processorAB2BoundaryIntegration::correctPStruct(real dt, const realx3PointField_D &vel)
+{
+
+	#ifndef BoundaryModel1
+    if(this->isBoundaryMaster())
+    {
+		const uint32 thisIndex = thisBoundaryIndex();
+		const auto&  AB2 =  static_cast<const AdamsBashforth2&>(Integration());
+		const auto& dy1View = AB2.BoundaryField(thisIndex).neighborProcField().deviceView();
+		const auto& velView = vel.BoundaryField(thisIndex).neighborProcField().deviceView();
+		const auto& xposView = boundary().neighborProcPoints().deviceView();        
+		const rangeU32 aRange(0u, dy1View.size());
+		return  AB2Kernels::intAllActive(
+			"AB2Integration::correctPStruct."+this->boundaryName(), 
+			dt,
+			aRange, 
+			xposView, 
+			velView, 
+			dy1View
+		); 
+    }
+#endif //BoundaryModel1
+    
+    return true;
+}

src/Integration/AdamsBashforth2/processorAB2BoundaryIntegration.hpp

@@ -0,0 +1,51 @@
+
+
+#ifndef __processorAB2BoundaryIntegration_hpp__
+#define __processorAB2BoundaryIntegration_hpp__
+
+#include "boundaryIntegration.hpp"
+
+namespace pFlow
+{
+
+class processorAB2BoundaryIntegration
+:
+    public boundaryIntegration
+{   
+public:
+
+    TypeInfo("boundaryIntegration<processor,AdamsBashforth2>");
+
+    processorAB2BoundaryIntegration(
+        const boundaryBase& boundary,
+        const pointStructure& pStruct, 
+        const word& method,
+        integration& intgrtn
+    );
+
+    ~processorAB2BoundaryIntegration()override=default;
+
+    
+	bool correct(
+        real dt, 
+        const realx3PointField_D& y, 
+        const realx3PointField_D& dy)override;
+	
+
+    
+    bool correctPStruct(real dt, const realx3PointField_D& vel)override;
+    
+    
+    add_vCtor(
+        boundaryIntegration,
+        processorAB2BoundaryIntegration,
+        boundaryBase
+    );
+
+    
+
+};
+
+}
+
+#endif

src/Interaction/contactSearch/boundaries/processorBoundaryContactSearch/processorBoundaryContactSearch.cpp

@@ -0,0 +1,111 @@
+/*------------------------------- 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 "processorBoundaryContactSearch.hpp"
+#include "contactSearch.hpp"
+#include "particles.hpp"
+//#include "pointStructure.hpp"
+//#include "geometry.hpp"
+
+
+void pFlow::processorBoundaryContactSearch::setSearchBox()
+{
+       
+    auto l = boundary().neighborLength();
+    auto n = boundary().boundaryPlane().normal();
+    auto pp1 = boundary().boundaryPlane().parallelPlane(l);
+    auto pp2 = boundary().boundaryPlane().parallelPlane(-l);
+
+    realx3 minP1 = min(min(min(pp1.p1(), pp1.p2()), pp1.p3()), pp1.p4());
+    realx3 maxP1 = max(max(max(pp1.p1(), pp1.p2()), pp1.p3()), pp1.p4());
+
+    realx3 minP2 = min(min(min(pp2.p1(), pp2.p2()), pp2.p3()), pp2.p4());
+    realx3 maxP2 = max(max(max(pp2.p1(), pp2.p2()), pp2.p3()), pp2.p4());
+
+    auto minP = min(minP1, minP2) - l*(realx3(1.0)-abs(n));
+    auto maxP = max(maxP1, maxP2) + l*(realx3(1.0)-abs(n));
+    
+    searchBox_={minP, maxP};
+}
+
+pFlow::processorBoundaryContactSearch::processorBoundaryContactSearch(
+    const dictionary &dict,
+    const boundaryBase &boundary,
+    const contactSearch &cSearch)
+: 
+    boundaryContactSearch(dict, boundary, cSearch),
+    diameter_(cSearch.Particles().boundingSphere()),
+    masterSearch_(this->isBoundaryMaster()),
+    sizeRatio_(dict.getVal<real>("sizeRatio"))
+{
+    
+    if(masterSearch_)
+    {        
+        setSearchBox();
+        
+        real minD;
+        real maxD;
+        cSearch.Particles().boundingSphereMinMax(minD, maxD);
+        
+        ppContactSearch_ = makeUnique<twoPartContactSearch>(
+            searchBox_,
+            maxD,
+            sizeRatio_);
+    }
+    else
+    {
+        searchBox_={{0,0,0},{0,0,0}};
+    }
+}
+
+bool pFlow::processorBoundaryContactSearch::broadSearch
+(
+    uint32 iter, 
+    real t, 
+    real dt, 
+    csPairContainerType &ppPairs, 
+    csPairContainerType &pwPairs, 
+    bool force
+)
+{
+    if(masterSearch_)
+    {
+        const auto thisPoints = boundary().thisPoints();
+        const auto& neighborProcPoints = boundary().neighborProcPoints();
+        const auto& bDiams = diameter_.BoundaryField(thisBoundaryIndex());
+        const auto thisDiams = bDiams.thisField();
+        const auto& neighborProcDiams = bDiams.neighborProcField();
+        
+        ppContactSearch_().broadSearchPP(
+            ppPairs, 
+            thisPoints, 
+            thisDiams,
+            neighborProcPoints,
+            neighborProcDiams,
+            boundaryName()
+        );
+        //pOutput<<"ppSize "<< ppPairs.size()<<endl;
+        return true;
+
+    }else
+    {
+        return true;
+    }
+}

src/Interaction/contactSearch/boundaries/processorBoundaryContactSearch/processorBoundaryContactSearch.hpp

@@ -0,0 +1,76 @@
+/*------------------------------- 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 __processorBoundaryContactSearch_hpp__
+#define __processorBoundaryContactSearch_hpp__
+
+#include "boundaryContactSearch.hpp"
+#include "pointFields.hpp"
+#include "twoPartContactSearch.hpp"
+
+namespace pFlow
+{
+
+class processorBoundaryContactSearch : public boundaryContactSearch
+{
+private:
+
+	box                             searchBox_;
+
+	uniquePtr<twoPartContactSearch> ppContactSearch_ = nullptr;
+
+	const realPointField_D&         diameter_;
+
+	bool                            masterSearch_;
+
+	real 							sizeRatio_;
+
+	void                            setSearchBox();
+
+public:
+
+	TypeInfo("boundaryContactSearch<MPI,processor>")
+
+	  processorBoundaryContactSearch(
+	    const dictionary&    dict,
+	    const boundaryBase&  boundary,
+	    const contactSearch& cSearch
+	  );
+
+	~processorBoundaryContactSearch() override = default;
+
+	add_vCtor(
+	  boundaryContactSearch,
+	  processorBoundaryContactSearch,
+	  boundaryBase
+	);
+
+	bool broadSearch(
+	  uint32               iter,
+	  real                 t,
+	  real                 dt,
+	  csPairContainerType& ppPairs,
+	  csPairContainerType& pwPairs,
+	  bool                 force = false
+	) override;
+};
+
+}
+
+#endif //__processorBoundaryContactSearch_hpp__

src/Interaction/contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp

@@ -0,0 +1,163 @@
+
+#include "twoPartContactSearch.hpp"
+#include "twoPartContactSearchKernels.hpp"
+#include "phasicFlowKokkos.hpp"
+#include "streams.hpp"
+
+void pFlow::twoPartContactSearch::checkAllocateNext(uint32 n)
+{
+	if( nextCapacity_ < n)
+	{
+		nextCapacity_ = n;	
+		reallocNoInit(next_, n);
+	}
+}
+
+void pFlow::twoPartContactSearch::nullifyHead()
+{
+	fill(head_, static_cast<uint32>(-1));
+}
+
+void pFlow::twoPartContactSearch::nullifyNext(uint32 n)
+{
+	fill(next_, 0u, n, static_cast<uint32>(-1));
+}
+
+void pFlow::twoPartContactSearch::buildList(
+	const deviceScatteredFieldAccess<realx3> &points)
+{
+	if(points.empty())return;
+	uint32 n = points.size();
+	checkAllocateNext(n);
+	nullifyNext(n);
+	nullifyHead();
+	
+	pFlow::twoPartContactSearchKernels::buildNextHead(
+		points,
+		searchCells_,
+		head_,
+		next_
+	);
+}
+
+pFlow::twoPartContactSearch::twoPartContactSearch
+(
+    const box &domain,
+    real cellSize,
+	real sizeRatio
+)
+:
+	searchCells_(domain, cellSize),
+	head_("periodic:head",searchCells_.nx(), searchCells_.ny(), searchCells_.nz()),
+	sizeRatio_(sizeRatio)
+{
+
+}
+
+bool pFlow::twoPartContactSearch::broadSearchPP
+(
+	csPairContainerType &ppPairs, 
+	const deviceScatteredFieldAccess<realx3> &points1, 
+	const deviceScatteredFieldAccess<real>& diams1, 
+	const deviceScatteredFieldAccess<realx3> &points2, 
+	const deviceScatteredFieldAccess<real>& diams2, 
+	const realx3& transferVec
+)
+{
+	if(points1.empty())return true;
+	if(points2.empty()) return true;
+
+	buildList(points1);
+
+	uint32 nNotInserted = 1;
+
+	// loop until the container size fits the numebr of contact pairs
+	while (nNotInserted > 0)
+	{
+
+		nNotInserted = pFlow::twoPartContactSearchKernels::broadSearchPP
+		(
+			ppPairs,
+			points1,
+			diams1,
+			points2,
+			diams2,
+			transferVec,
+			head_,
+			next_,
+			searchCells_,
+			sizeRatio_
+		);
+		
+	
+		if(nNotInserted)
+		{
+			// - resize the container
+			//   note that getFull now shows the number of failed insertions.
+			uint32 len = max(nNotInserted,100u) ;
+			
+			auto oldCap = ppPairs.capacity();
+			
+			ppPairs.increaseCapacityBy(len);
+
+			INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
+			oldCap << " to "<<ppPairs.capacity()<<" in contact search in boundary region."<<END_INFO;
+			
+		}
+
+	}	
+	
+    return true;
+}
+
+bool pFlow::twoPartContactSearch::broadSearchPP
+(
+	csPairContainerType &ppPairs,
+	const deviceScatteredFieldAccess<realx3> &points1,
+	const deviceScatteredFieldAccess<real> &diams1,
+	const realx3Vector_D& points2,
+	const realVector_D& diams2,
+	const word& name
+)
+{
+	buildList(points1);
+
+	uint32 nNotInserted = 1;
+
+	// loop until the container size fits the numebr of contact pairs
+	while (nNotInserted > 0)
+	{
+
+		nNotInserted = pFlow::twoPartContactSearchKernels::broadSearchPP
+		(
+			ppPairs,
+			points1,
+			diams1,
+			points2,
+			diams2,
+			head_,
+			next_,
+			searchCells_,
+			sizeRatio_
+		);
+		
+	
+		if(nNotInserted)
+		{
+			// - resize the container
+			//   note that getFull now shows the number of failed insertions.
+			uint32 len = max(nNotInserted,100u) ;
+			
+			auto oldCap = ppPairs.capacity();
+			
+			ppPairs.increaseCapacityBy(len);
+			
+			INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
+			oldCap << " to "<<ppPairs.capacity()<<" in boundary contact search in "<< name <<END_INFO;
+			
+		}
+
+	}	
+	
+    return true;
+}

src/Interaction/contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.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 __twoPartContactSearch_hpp__
+#define __twoPartContactSearch_hpp__
+
+#include "contactSearchGlobals.hpp"
+#include "scatteredFieldAccess.hpp"
+#include "cells.hpp"
+#include "VectorSingles.hpp"
+
+namespace pFlow
+{
+
+class twoPartContactSearch
+{
+public:
+	using HeadType = deviceViewType3D<uint32>;
+
+	using NextType = deviceViewType1D<uint32>;
+
+private:
+
+	cells    searchCells_;
+
+	HeadType head_{ "periodic::head", 1, 1, 1 };
+
+	NextType next_{ "periodic::next", 1 };
+
+	real     sizeRatio_ = 1.0;
+
+	uint32   nextCapacity_ = 0;
+
+	void checkAllocateNext(uint32 n);
+
+	void nullifyHead();
+
+	void nullifyNext(uint32 n);
+
+	void buildList(
+		const deviceScatteredFieldAccess<realx3> &points);
+
+public:
+	twoPartContactSearch(
+		const box &domain,
+		real cellSize,
+		real sizeRatio = 1.0);
+
+	/// @brief Perform a broad-search for spheres in two adjacent regions.
+	/// Region 1 is considered as the master (primary) region and region 2 as slave 
+	/// @param ppPairs pairs container which holds i and j
+	/// @param points1 point positions in region 1
+	/// @param diams1  diameter of spheres in region 1
+	/// @param points2 point positions in region 2
+	/// @param diams2 diameter of spheres in region 2
+	/// @param transferVec a vector to transfer points from region 2 to region 1
+	/// @return true if it is successful 
+	bool broadSearchPP(
+		csPairContainerType &ppPairs,
+		const deviceScatteredFieldAccess<realx3> &points1,
+		const deviceScatteredFieldAccess<real> &diams1,
+		const deviceScatteredFieldAccess<realx3> &points2,
+		const deviceScatteredFieldAccess<real> &diams2,
+		const realx3 &transferVec);
+	
+	bool broadSearchPP(
+		csPairContainerType &ppPairs,
+		const deviceScatteredFieldAccess<realx3> &points1,
+		const deviceScatteredFieldAccess<real> &diams1,
+		const realx3Vector_D& points2,
+		const realVector_D& diams2,
+		const word& name);
+
+	const auto& searchCells()const
+	{
+		return searchCells_;
+	}
+
+	real sizeRatio()const
+	{
+		return sizeRatio_;
+	}
+};
+
+}
+
+#endif //__twoPartContactSearch_hpp__

src/Interaction/contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp

@@ -0,0 +1,186 @@
+#include "twoPartContactSearchKernels.hpp"
+
+INLINE_FUNCTION_HD
+bool
+sphereSphereCheckB(
+  const pFlow::realx3& p1,
+  const pFlow::realx3  p2,
+  pFlow::real          d1,
+  pFlow::real          d2
+)
+{
+	return pFlow::length(p2 - p1) < 0.5 * (d2 + d1);
+}
+
+void
+pFlow::twoPartContactSearchKernels::buildNextHead(
+  const deviceScatteredFieldAccess<realx3>& points,
+  const cells&                              searchCells,
+  deviceViewType3D<uint32>&                 head,
+  deviceViewType1D<uint32>&                 next
+)
+{
+	
+	uint32 n = points.size();
+
+	Kokkos::parallel_for(
+	  "pFlow::ppwBndryContactSearch::buildList",
+	  deviceRPolicyStatic(0, n),
+	  LAMBDA_HD(uint32 i) {
+		  int32x3 ind;
+		  if (searchCells.pointIndexInDomain(points[i], ind))
+		  {
+			  // discards points out of searchCell
+			  uint32 old =
+			    Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
+			  next[i] = old;
+		  }
+	  }
+	);
+	Kokkos::fence();
+}
+
+pFlow::uint32
+pFlow::twoPartContactSearchKernels::broadSearchPP(
+  csPairContainerType&                      ppPairs,
+  const deviceScatteredFieldAccess<realx3>& points,
+  const deviceScatteredFieldAccess<real>&   diams,
+  const deviceScatteredFieldAccess<realx3>& mirrorPoints,
+  const deviceScatteredFieldAccess<real>&   mirrorDiams,
+  const realx3&                             transferVec,
+  const deviceViewType3D<uint32>&           head,
+  const deviceViewType1D<uint32>&           next,
+  const cells&                              searchCells,
+  const real                                sizeRatio
+)
+{
+	if (points.empty())
+		return 0;
+	if (mirrorPoints.empty())
+		return 0;
+
+	auto   nMirror = mirrorPoints.size();
+
+	uint32 getFull = 0;
+
+	Kokkos::parallel_reduce(
+	  "pFlow::twoPartContactSearchKernels::broadSearchPP",
+	  deviceRPolicyStatic(0, nMirror),
+	  LAMBDA_HD(const uint32 mrrI, uint32& getFullUpdate) {
+		  realx3  p_m = mirrorPoints(mrrI) + transferVec;
+
+		  int32x3 ind_m;
+		  if (!searchCells.pointIndexInDomain(p_m, ind_m))
+			  return;
+
+		  real d_m = sizeRatio * mirrorDiams[mrrI];
+
+		  for (int ii = -1; ii < 2; ii++)
+		  {
+			  for (int jj = -1; jj < 2; jj++)
+			  {
+				  for (int kk = -1; kk < 2; kk++)
+				  {
+					  auto ind = ind_m + int32x3{ ii, jj, kk };
+
+					  if (!searchCells.inCellRange(ind))
+						  continue;
+
+					  uint32 thisI = head(ind.x(), ind.y(), ind.z());
+					  while (thisI !=  static_cast<uint32>(-1))
+					  {
+						  auto d_n = sizeRatio * diams[thisI];
+
+						  // first item is for this boundary and second itme,
+						  // for mirror
+						  if(sphereSphereCheckB(p_m, points[thisI], d_m, d_n)&&
+                       ppPairs.insert(thisI,mrrI) == static_cast<uint32>(-1))
+						  {
+							  getFullUpdate++;
+						  }
+
+						  thisI = next(thisI);
+					  }
+				  }
+			  }
+		  }
+	  },
+	  getFull
+	);
+
+	return getFull;
+}
+
+pFlow::uint32
+pFlow::twoPartContactSearchKernels::broadSearchPP(
+  csPairContainerType&                      ppPairs,
+  const deviceScatteredFieldAccess<realx3>& points1,
+  const deviceScatteredFieldAccess<real>&   diams1,
+  const realx3Vector_D&                     points2,
+  const realVector_D&                       diams2,
+  const deviceViewType3D<uint32>&           head,
+  const deviceViewType1D<uint32>&           next,
+  const cells&                              searchCells,
+  real                                      sizeRatio
+)
+{
+	if (points1.empty())
+		return 0;
+	if (points2.empty())
+		return 0;
+
+	auto   nP2         = points2.size();
+	auto   points2View = points2.deviceView();
+	auto   diams2View  = diams2.deviceView();
+
+	uint32 getFull = 0;
+
+	Kokkos::parallel_reduce(
+	  "pFlow::twoPartContactSearchKernels::broadSearchPP",
+	  deviceRPolicyStatic(0, nP2),
+	  LAMBDA_HD(const uint32 i2, uint32& getFullUpdate) {
+		  realx3  p_m = points2View(i2);
+
+		  int32x3 ind_m;
+		  if (!searchCells.pointIndexInDomain(p_m, ind_m))
+			  return;
+
+		  real d_m = sizeRatio * diams2View[i2];
+
+		  for (int ii = -1; ii < 2; ii++)
+		  {
+			  for (int jj = -1; jj < 2; jj++)
+			  {
+				  for (int kk = -1; kk < 2; kk++)
+				  {
+					  auto ind = ind_m + int32x3{ ii, jj, kk };
+
+					  if (!searchCells.inCellRange(ind))
+					  {
+						  continue;
+					  }
+
+					  uint32 i1 = head(ind.x(), ind.y(), ind.z());
+					  while (i1 != static_cast<uint32>(-1))
+					  {
+						  auto d_n = sizeRatio * diams1[i1];
+
+						  // first item is for this boundary and second itme,
+						  // for mirror
+						  if(sphereSphereCheckB(p_m, points1[i1], d_m, d_n)&&
+                       ppPairs.insert(i1,i2) == static_cast<uint32>(-1))
+						  {
+							  getFullUpdate++;
+						  }
+
+						  i1 = next(i1);
+					  }
+				  }
+			  }
+		  }
+	  },
+	  getFull
+	);
+
+	return getFull;
+}

src/Interaction/contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.hpp

@@ -0,0 +1,49 @@
+#ifndef __twoPartContactSearchKernels_hpp__
+#define __twoPartContactSearchKernels_hpp__
+
+#include "contactSearchGlobals.hpp"
+#include "cells.hpp"
+#include "contactSearchFunctions.hpp"
+#include "scatteredFieldAccess.hpp"
+#include "VectorSingles.hpp"
+
+namespace pFlow::twoPartContactSearchKernels
+{
+
+void buildNextHead(
+    const deviceScatteredFieldAccess<realx3> &points,
+    const cells              &searchCells,
+    deviceViewType3D<uint32> &head, 
+    deviceViewType1D<uint32> &next );
+
+
+uint32 broadSearchPP
+(
+	csPairContainerType 						&ppPairs,
+	const deviceScatteredFieldAccess<realx3> 	&points, 
+	const deviceScatteredFieldAccess<real> 		&diams, 
+	const deviceScatteredFieldAccess<realx3>	&mirrorPoints, 
+	const deviceScatteredFieldAccess<real> 		&mirrorDiams, 
+	const realx3 								&transferVec,
+    const deviceViewType3D<uint32>				&head,
+    const deviceViewType1D<uint32>				&next,
+    const cells									&searchCells,
+	real sizeRatio
+);
+
+uint32
+broadSearchPP(
+	csPairContainerType&                      ppPairs,
+	const deviceScatteredFieldAccess<realx3>& points1,
+	const deviceScatteredFieldAccess<real>&   diams1,
+	const realx3Vector_D&                     points2,
+	const realVector_D&                       diams2,
+	const deviceViewType3D<uint32>&           head,
+	const deviceViewType1D<uint32>&           next,
+	const cells&                              searchCells,
+	real                                      sizeRatio
+);
+}
+
+
+#endif //__twoPartContactSearchKernels_hpp__

src/Interaction/sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySIKernels.hpp

@@ -0,0 +1,132 @@
+
+#ifndef __processorBoundarySIKernels_hpp__
+#define __processorBoundarySIKernels_hpp__
+
+namespace pFlow::MPI::processorBoundarySIKernels
+{
+
+template<typename ContactListType, typename ContactForceModel>
+inline
+void sphereSphereInteraction
+(
+    const word& kernalName,
+    real dt,
+    const ContactListType&              cntctList,
+    const ContactForceModel&            forceModel,
+    const deviceScatteredFieldAccess<realx3>& thisPoints,
+    const deviceViewType1D<real>& thisDiam,
+    const deviceViewType1D<uint32>& thisPropId,
+    const deviceViewType1D<realx3>& thisVel,
+    const deviceViewType1D<realx3>& thisRVel,
+    const deviceViewType1D<realx3>& thisCForce,
+    const deviceViewType1D<realx3>& thisCTorque,
+    const deviceViewType1D<realx3>& neighborPoints,
+    const deviceViewType1D<real>& neighborDiam,
+    const deviceViewType1D<uint32>& neighborPropId,
+    const deviceViewType1D<realx3>& neighborVel,
+    const deviceViewType1D<realx3>& neighborRVel,
+    const deviceViewType1D<realx3>& neighborCForce,
+    const deviceViewType1D<realx3>& neighborCTorque
+)
+{
+    
+    using ValueType = typename ContactListType::ValueType;
+    uint32 ss = cntctList.size();
+    if(ss == 0u)return;
+
+    uint32 lastItem = cntctList.loopCount();
+
+    Kokkos::parallel_for(
+        kernalName,
+        deviceRPolicyDynamic(0,lastItem),
+        LAMBDA_HD(uint32 n)
+        {
+
+        if(!cntctList.isValid(n))return;
+
+            auto [i,j] = cntctList.getPair(n);
+            uint32 ind_i = thisPoints.index(i);
+            uint32 ind_j = j;
+
+            real Ri = 0.5*thisDiam[ind_i];
+            real Rj = 0.5*neighborDiam[ind_j];
+            realx3 xi = thisPoints.field()[ind_i];
+            realx3 xj = neighborPoints[ind_j];
+
+            real dist = length(xj-xi);
+            real ovrlp = (Ri+Rj) - dist;
+            
+            if( ovrlp >0.0 )
+            {
+                auto Nij = (xj-xi)/max(dist,smallValue);
+                auto wi = thisRVel[ind_i];
+                auto wj = neighborRVel[ind_j];
+                auto Vr = thisVel[ind_i] - neighborVel[ind_j] + cross((Ri*wi+Rj*wj), Nij);
+                
+                auto history = cntctList.getValue(n);
+
+                int32 propId_i = thisPropId[ind_i];
+                int32 propId_j = neighborPropId[ind_j];
+
+                realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
+
+                // calculates contact force 
+                forceModel.contactForce(
+                    dt, i, j,
+                    propId_i, propId_j,
+                    Ri, Rj,
+                    ovrlp,
+                    Vr, Nij,
+                    history,
+                    FCn, FCt);
+
+                forceModel.rollingFriction(
+                    dt, i, j,
+                    propId_i, propId_j,
+                    Ri, Rj,
+                    wi, wj,
+                    Nij,
+                    FCn,
+                    Mri, Mrj);
+
+                auto M = cross(Nij,FCt);
+                Mij = Ri*M+Mri;
+                Mji = Rj*M+Mrj;
+                
+                auto FC = FCn + FCt;
+                
+
+                Kokkos::atomic_add(&thisCForce[ind_i].x_,FC.x_);
+                Kokkos::atomic_add(&thisCForce[ind_i].y_,FC.y_);
+                Kokkos::atomic_add(&thisCForce[ind_i].z_,FC.z_);
+
+                Kokkos::atomic_add(&neighborCForce[ind_j].x_,-FC.x_);
+                Kokkos::atomic_add(&neighborCForce[ind_j].y_,-FC.y_);
+                Kokkos::atomic_add(&neighborCForce[ind_j].z_,-FC.z_);
+
+                Kokkos::atomic_add(&thisCTorque[ind_i].x_, Mij.x_);
+                Kokkos::atomic_add(&thisCTorque[ind_i].y_, Mij.y_);
+                Kokkos::atomic_add(&thisCTorque[ind_i].z_, Mij.z_);
+
+                Kokkos::atomic_add(&neighborCTorque[ind_j].x_, Mji.x_);
+                Kokkos::atomic_add(&neighborCTorque[ind_j].y_, Mji.y_);
+                Kokkos::atomic_add(&neighborCTorque[ind_j].z_, Mji.z_);
+                
+
+                cntctList.setValue(n,history);
+
+            }
+            else
+            {
+                cntctList.setValue(n, ValueType());
+            }
+
+        });
+    Kokkos::fence();
+}
+
+
+} //pFlow::MPI::processorBoundarySIKernels
+
+
+#endif //__processorBoundarySIKernels_hpp__

src/Interaction/sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteraction.cpp

@@ -0,0 +1,256 @@
+/*------------------------------- 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 "processorBoundarySIKernels.hpp"
+
+template <typename cFM, typename gMM>
+pFlow::MPI::processorBoundarySphereInteraction<cFM, gMM>::processorBoundarySphereInteraction(
+	const boundaryBase &boundary,
+	const sphereParticles &sphPrtcls,
+	const GeometryMotionModel &geomMotion)
+: 
+	boundarySphereInteraction<cFM,gMM>(
+		boundary, 
+		sphPrtcls, 
+		geomMotion
+	),
+	masterInteraction_(boundary.isBoundaryMaster())
+{
+	if(masterInteraction_)
+	{
+		this->allocatePPPairs();
+		this->allocatePWPairs();
+	}
+		
+}
+
+
+#ifdef BoundaryModel1
+
+template <typename cFM, typename gMM>
+bool pFlow::MPI::processorBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
+(
+	real dt,
+	const ContactForceModel &cfModel,
+	uint32 step
+)
+{
+	
+	// master processor calculates the contact force/torque and sends data back to the 
+	// neighbor processor (slave processor).
+	// slave processor recieves the data and adds the data to the internalField 
+	if(masterInteraction_) 
+	{
+		if(step==1)return true;
+
+		const auto & sphPar = this->sphParticles();
+		uint32 thisIndex = this->boundary().thisBoundaryIndex();
+		
+		const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&> (
+		sphPar.contactForce().BoundaryField(thisIndex));
+
+		const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+		sphPar.contactTorque().BoundaryField(thisIndex));
+
+		if(step == 2 )
+		{
+			iter++;
+			pFlow::MPI::processorBoundarySIKernels::sphereSphereInteraction(
+				"ppBoundaryInteraction."+this->boundaryName(),
+				dt,
+				this->ppPairs(),
+				cfModel,
+				this->boundary().thisPoints(),
+				sphPar.diameter().deviceViewAll(),
+				sphPar.propertyId().deviceViewAll(),
+				sphPar.velocity().deviceViewAll(),
+				sphPar.rVelocity().deviceViewAll(),
+				sphPar.contactForce().deviceViewAll(),
+				sphPar.contactTorque().deviceViewAll(),
+				this->boundary().neighborProcPoints().deviceViewAll(),
+				sphPar.diameter().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.propertyId().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.velocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.rVelocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				cfBndry.neighborProcField().deviceViewAll(),
+				ctBndry.neighborProcField().deviceViewAll()
+			);
+			
+			return true;
+		}		
+		else if(step == 3 )
+		{
+			cfBndry.sendBackData();
+			ctBndry.sendBackData();
+			
+			return true;
+		}
+		
+		
+		return false;	
+	}
+	else
+	{
+		
+		if(step == 1 )
+		{
+			const auto & sphPar = this->sphParticles();
+			uint32 thisIndex = this->boundary().thisBoundaryIndex();
+			const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&>(
+			sphPar.contactForce().BoundaryField(thisIndex));
+			const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+			sphPar.contactTorque().BoundaryField(thisIndex));
+
+			cfBndry.recieveBackData();
+			ctBndry.recieveBackData();
+			
+			return false;
+		}
+		else if(step == 11)
+		{
+			const auto & sphPar = this->sphParticles();
+			uint32 thisIndex = this->boundary().thisBoundaryIndex();
+			const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&>(
+			sphPar.contactForce().BoundaryField(thisIndex));
+			const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+			sphPar.contactTorque().BoundaryField(thisIndex));
+
+			cfBndry.addBufferToInternalField();
+			ctBndry.addBufferToInternalField();
+	
+			return true;
+		}
+		
+		
+		return false;
+	}
+
+	return false;
+}
+#else
+
+template <typename cFM, typename gMM>
+bool pFlow::MPI::processorBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
+(
+	real dt,
+	const ContactForceModel &cfModel,
+	uint32 step
+)
+{
+	
+	// master processor calculates the contact force/torque and sends data back to the 
+	// neighbor processor (slave processor).
+	// slave processor recieves the data and adds the data to the internalField 
+	if(masterInteraction_) 
+	{
+		if(step==1)return true;
+
+		const auto & sphPar = this->sphParticles();
+		uint32 thisIndex = this->boundary().thisBoundaryIndex();
+		
+		const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&> (
+		sphPar.contactForce().BoundaryField(thisIndex));
+
+		const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+		sphPar.contactTorque().BoundaryField(thisIndex));
+
+		if(step == 2 )
+		{
+			
+			pFlow::MPI::processorBoundarySIKernels::sphereSphereInteraction(
+				"ppBoundaryInteraction."+this->boundaryName(),
+				dt,
+				this->ppPairs(),
+				cfModel,
+				this->boundary().thisPoints(),
+				sphPar.diameter().deviceViewAll(),
+				sphPar.propertyId().deviceViewAll(),
+				sphPar.velocity().deviceViewAll(),
+				sphPar.rVelocity().deviceViewAll(),
+				sphPar.contactForce().deviceViewAll(),
+				sphPar.contactTorque().deviceViewAll(),
+				this->boundary().neighborProcPoints().deviceViewAll(),
+				sphPar.diameter().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.propertyId().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.velocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				sphPar.rVelocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
+				cfBndry.neighborProcField().deviceViewAll(),
+				ctBndry.neighborProcField().deviceViewAll()
+			);
+		
+			return true;
+		}		
+		else if(step == 3 )
+		{
+			cfBndry.sendBackData();
+			ctBndry.sendBackData();
+			return true;
+		}
+		else if(step == 11 )
+		{
+			cfBndry.updateBoundaryFromSlave();
+			ctBndry.updateBoundaryFromSlave();
+			return true;
+		}
+
+		return false;	
+	}
+	else
+	{
+		
+		if(step == 1 )
+		{
+			const auto & sphPar = this->sphParticles();
+			uint32 thisIndex = this->boundary().thisBoundaryIndex();
+			const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&>(
+			sphPar.contactForce().BoundaryField(thisIndex));
+			const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+			sphPar.contactTorque().BoundaryField(thisIndex));
+
+			cfBndry.recieveBackData();
+			ctBndry.recieveBackData();
+			
+			return false;
+		}
+		else if(step == 11)
+		{
+			const auto & sphPar = this->sphParticles();
+			uint32 thisIndex = this->boundary().thisBoundaryIndex();
+			const auto& cfBndry = static_cast<const processorBoundaryField<realx3>&>(
+			sphPar.contactForce().BoundaryField(thisIndex));
+			const auto& ctBndry = static_cast<const processorBoundaryField<realx3>&> (
+			sphPar.contactTorque().BoundaryField(thisIndex));
+
+			cfBndry.addBufferToInternalField();
+			cfBndry.updateBoundaryToMaster();
+
+			ctBndry.addBufferToInternalField();
+			ctBndry.updateBoundaryToMaster();
+
+			return true;
+		}
+				
+		return false;
+	}
+
+	return false;
+}
+
+#endif

src/Interaction/sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteraction.hpp

@@ -0,0 +1,93 @@
+/*------------------------------- 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 __processorBoundarySphereInteraction_hpp__
+#define __processorBoundarySphereInteraction_hpp__
+
+#include "boundarySphereInteraction.hpp"
+#include "processorBoundaryField.hpp"
+#include "boundaryProcessor.hpp"
+
+namespace pFlow::MPI
+{
+
+template<typename contactForceModel,typename geometryMotionModel>
+class processorBoundarySphereInteraction
+:
+    public boundarySphereInteraction<contactForceModel, geometryMotionModel>
+{
+public:
+
+    using PBSInteractionType = 
+        processorBoundarySphereInteraction<contactForceModel,geometryMotionModel>;
+
+    using BSInteractionType = 
+        boundarySphereInteraction<contactForceModel, geometryMotionModel>;
+
+    using GeometryMotionModel 	= typename BSInteractionType::GeometryMotionModel;
+	
+	using ContactForceModel 	= typename BSInteractionType::ContactForceModel;
+
+	using MotionModel 			= typename geometryMotionModel::MotionModel;
+
+	using ModelStorage 			= typename ContactForceModel::contactForceStorage;
+	
+	using IdType 			    = typename BSInteractionType::IdType;
+
+	using IndexType    		    = typename BSInteractionType::IndexType;
+    
+    using ContactListType       = typename BSInteractionType::ContactListType; 
+
+private:
+    
+    bool    masterInteraction_;
+
+public:
+
+    TypeInfoTemplate22("boundarySphereInteraction", "processor",ContactForceModel, MotionModel);
+
+	
+    processorBoundarySphereInteraction(
+        const boundaryBase& boundary,
+		const sphereParticles& sphPrtcls,
+		const GeometryMotionModel& geomMotion
+    );
+
+    add_vCtor
+    (
+        BSInteractionType,
+        PBSInteractionType,
+        boundaryBase
+    );
+    
+    ~processorBoundarySphereInteraction()override = default;
+
+    bool sphereSphereInteraction(
+        real dt,
+		const ContactForceModel& cfModel,
+        uint32 step)override;
+	
+};
+
+}
+
+#include "processorBoundarySphereInteraction.cpp"
+
+
+#endif //__processorBoundarySphereInteraction_hpp__

src/Interaction/sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteractions.cpp

@@ -0,0 +1,17 @@
+
+#include "processorBoundarySphereInteraction.hpp"
+#include "geometryMotions.hpp"
+#include "contactForceModels.hpp"
+
+
+template class pFlow::MPI::processorBoundarySphereInteraction
+<
+    pFlow::cfModels::limitedNonLinearNormalRolling,
+    pFlow::rotationAxisMotionGeometry
+>;
+
+template class pFlow::MPI::processorBoundarySphereInteraction
+<
+    pFlow::cfModels::nonLimitedNonLinearNormalRolling,
+    pFlow::rotationAxisMotionGeometry
+>;

src/Interaction/sphereInteraction/sphereInteraction/sphereInteraction.cpp

@@ -359,7 +359,7 @@ bool pFlow::sphereInteraction<cFM,gMM, cLT>::hearChanges
 	if(msg.equivalentTo(message::ITEMS_REARRANGE))
 	{
 		notImplementedFunction;
-		return false;
+		return true;
 	}
 	
 	fatalErrorInFunction<<"Event "<< msg.eventNames()<<

src/Particles/SphereParticles/processorBoundarySphereParticles.cpp

@@ -0,0 +1,46 @@
+#include "processorBoundarySphereParticles.hpp"
+#include "sphereParticles.hpp"
+#include "boundaryProcessor.hpp"
+
+pFlow::processorBoundarySphereParticles::processorBoundarySphereParticles(
+    const boundaryBase &boundary, 
+    sphereParticles &prtcls
+)
+:
+    boundarySphereParticles(boundary, prtcls)
+{
+
+}
+
+bool pFlow::processorBoundarySphereParticles::acceleration(const timeInfo &ti, const realx3& g)
+{
+    
+
+#ifndef BoundaryModel1
+    
+
+    if(isBoundaryMaster())
+    {
+        auto thisIndex = thisBoundaryIndex();
+        auto mass = Particles().mass().BoundaryField(thisIndex).neighborProcField().deviceView();
+        auto I = Particles().I().BoundaryField(thisIndex).neighborProcField().deviceView();
+        auto cf = Particles().contactForce().BoundaryField(thisIndex).neighborProcField().deviceView();
+        auto ct = Particles().contactTorque().BoundaryField(thisIndex).neighborProcField().deviceView();
+        auto acc = Particles().acceleration().BoundaryField(thisIndex).neighborProcField().deviceView();
+        auto rAcc = Particles().rAcceleration().BoundaryField(thisIndex).neighborProcField().deviceView();
+
+        Kokkos::parallel_for(
+            "processorBoundary::acceleration."+this->boundaryName(),
+            deviceRPolicyStatic(0,mass.size()),
+            LAMBDA_HD(uint32 i){
+                acc[i] = cf[i]/mass[i] + g;
+                rAcc[i] = ct[i]/I[i];
+            });
+        Kokkos::fence();    
+    }
+
+#endif
+    
+
+    return true;
+}

src/Particles/SphereParticles/processorBoundarySphereParticles.hpp

@@ -0,0 +1,38 @@
+#ifndef __processorBoundarySphereParticles_hpp__
+#define __processorBoundarySphereParticles_hpp__
+
+#include "boundarySphereParticles.hpp"
+
+namespace pFlow
+{
+
+class processorBoundarySphereParticles
+:
+    public boundarySphereParticles
+{
+
+public:
+
+    /// type info
+    TypeInfo("boundarySphereParticles<MPI,processor>");
+
+    processorBoundarySphereParticles(
+        const boundaryBase &boundary,
+        sphereParticles&  prtcls
+    );
+
+    add_vCtor(
+        boundarySphereParticles,
+        processorBoundarySphereParticles,
+        boundaryBase
+    );
+
+    bool acceleration(const timeInfo& ti, const realx3& g)override;
+    
+
+};
+
+}
+
+
+#endif

src/Particles/particles/MPIParticleIdHandler/MPIParticleIdHandler.cpp

@@ -0,0 +1,70 @@
+#include "MPIParticleIdHandler.hpp"
+#include "procCommunication.hpp"
+
+pFlow::MPI::MPIParticleIdHandler::MPIParticleIdHandler
+(
+    pointStructure& pStruct
+)
+:
+    particleIdHandler(pStruct)
+{
+    initialIdCheck();
+}
+
+pFlow::Pair<pFlow::uint32, pFlow::uint32> 
+    pFlow::MPI::MPIParticleIdHandler::getIdRange(uint32 nNewParticles)
+{
+    uint32 startId;
+    if(maxId_==-1)
+    {
+        startId = 0;
+    }
+    else
+    {
+        startId = maxId_+1;
+    }
+    uint32 endId = startId+nNewParticles-1;
+    maxId_ = endId;
+    return {startId, endId};
+}
+
+bool pFlow::MPI::MPIParticleIdHandler::initialIdCheck()
+{
+    /// empty point structure / no particles in simulation
+    uint32 maxId = -1;
+    if( !pStruct().empty() ) 
+    {
+        maxId = max( *this );
+    }    
+    
+    auto maxIdAll = procVector<uint32>(pFlowProcessors());
+    auto numAll = procVector<uint32>(pFlowProcessors());
+    auto comm = procCommunication(pFlowProcessors());
+
+    comm.collectAllToAll(maxId, maxIdAll);
+    comm.collectAllToAll(size(),numAll);
+    
+    uint32 n = 0;
+    for(uint32 i=0; i<maxIdAll.size(); i++)
+    {
+        if( maxIdAll[i]==-1 && numAll[i]!= 0)
+        {
+            if(comm.localRank() == i)
+            {
+                fillSequence(*this, n);
+                maxId_ = size()-1 + n;
+            }
+            
+        }
+        else
+        {
+            if(comm.localRank() == i)
+            {
+                maxId_ = maxIdAll[i];
+            }
+        }
+        n += numAll[i];
+    }
+    
+    return true;
+}

src/Particles/particles/MPIParticleIdHandler/MPIParticleIdHandler.hpp

@@ -0,0 +1,60 @@
+/*------------------------------- 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 __MPIParticleIdHandler_hpp__
+#define __MPIParticleIdHandler_hpp__
+
+#include "particleIdHandler.hpp"
+
+namespace pFlow::MPI
+{
+
+class MPIParticleIdHandler : public particleIdHandler
+{
+private:
+
+	uint32 maxId_ = -1;
+
+	bool   initialIdCheck() override;
+
+public:
+
+	ClassInfo("particleIdHandler<MPI>");
+
+	explicit MPIParticleIdHandler(pointStructure& pStruct);
+
+	~MPIParticleIdHandler() override = default;
+
+	add_vCtor(
+    particleIdHandler, 
+    MPIParticleIdHandler, 
+    pointStructure
+  );
+
+	Pair<uint32, uint32> getIdRange(uint32 nNewParticles) override;
+
+	uint32               maxId() const override
+	{
+		return maxId_;
+	}
+};
+
+}
+
+#endif //__MPIParticleIdHandler_hpp__

src/Particles/particles/particles.hpp

@@ -185,6 +185,18 @@ public:
 		return contactTorque_;
 	}
 
+	inline 
+	uint32PointField_D& particleId()
+	{
+		return idHandler_();
+	}
+
+	inline
+	const uint32PointField_D& particleId() const
+	{
+		return idHandler_();
+	}
+
 	inline 
 	uint32 maxId()const
 	{

src/PostprocessData/CMakeLists.txt

@@ -9,9 +9,11 @@ set(SourceFiles
     # Regions
     region/regionPoints/regionPoints/regionPoints.cpp
     region/regionPoints/sphereRegionPoints/sphereRegionPoints.cpp
+    region/regionPoints/boxRegionPoints/boxRegionPoints.cpp
     region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
     region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
     region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
+    region/regionPoints/rectMeshRegionPoints/rectMeshRegionPoints.cpp
 
     # Postprocess components
     postprocessComponent/postprocessComponent/postprocessComponent.cpp

src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp

@@ -467,7 +467,7 @@ pFlow::postprocessData::fieldsDataBase::fieldsDataBase
     systemControl& control,
     const dictionary& postDict,
     bool inSimulation,
-    timeValue startTime
+    TimeValueType startTime
 )
 :
     time_(control.time()),
@@ -492,7 +492,7 @@ pFlow::postprocessData::fieldsDataBase::fieldsDataBase
     }
 }
 
-pFlow::timeValue pFlow::postprocessData::fieldsDataBase::currentTime() const
+pFlow::TimeValueType pFlow::postprocessData::fieldsDataBase::currentTime() const
 {
     return time_.currentTime();
 }
@@ -914,7 +914,7 @@ pFlow::uniquePtr<pFlow::postprocessData::fieldsDataBase>
     systemControl& control, 
     const dictionary& postDict, 
     bool inSimulation,
-    timeValue startTime
+    TimeValueType startTime
 )
 {
     word dbType;

src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp

@@ -78,7 +78,7 @@ private:
         anyList             allFields_;
 
         /// Map to store the last capture time of each field
-        wordMap<timeValue>  captureTime_;
+        wordMap<TimeValueType>  captureTime_;
 
         /// Reference to the Time object
         Time&               time_; 
@@ -178,7 +178,7 @@ public:
             systemControl& control, 
             const dictionary& postDict,
             bool inSimulation,
-            timeValue startTime);
+            TimeValueType startTime);
 
         /// no copy constructor
         fieldsDataBase(const fieldsDataBase&) = delete;
@@ -203,7 +203,7 @@ public:
                 systemControl& control, 
                 const dictionary& postDict, 
                 bool inSimulation,
-                timeValue startTime
+                TimeValueType startTime
             ),
             (control, postDict, inSimulation, startTime)
         );
@@ -211,7 +211,7 @@ public:
     
     // - Public Access Functions
         /// returns the current time
-        timeValue currentTime()const;
+        TimeValueType currentTime()const;
 
         /// const ref to object Time 
         const Time& time()const 
@@ -282,7 +282,7 @@ public:
         /// Get the next avaiable time folder after the current time folder 
         /// This is only used for post-simulation processing 
         virtual
-        timeValue getNextTimeFolder()const
+        TimeValueType getNextTimeFolder()const
         {
             return -1.0;
         }
@@ -291,7 +291,7 @@ public:
         /// This is used only for post-simulation processing
         /// @returns the time value of the next folder.
         virtual 
-        timeValue setToNextTimeFolder()
+        TimeValueType setToNextTimeFolder()
         {
             return -1.0;
         }
@@ -300,7 +300,7 @@ public:
         /// This is used only for post-simulation processing 
         /// @returns the time value of the skipped folder 
         virtual 
-        timeValue skipNextTimeFolder()
+        TimeValueType skipNextTimeFolder()
         {
             return -1.0;
         }
@@ -316,7 +316,7 @@ public:
         systemControl& control, 
         const dictionary& postDict, 
         bool inSimulation,
-        timeValue startTime);
+        TimeValueType startTime);
 };
 
 } // namespace pFlow::postprocessData

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.cpp

@@ -49,7 +49,7 @@ pFlow::postprocessData::simulationFieldsDataBase::simulationFieldsDataBase
     systemControl &control,
     const dictionary& postDict, 
     bool inSimulation,
-    timeValue startTime
+    TimeValueType startTime
 )
 :
     fieldsDataBase(control, postDict, inSimulation, startTime),

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.hpp

@@ -60,7 +60,7 @@ public:
         systemControl& control, 
         const dictionary& postDict,
         bool inSimulation,
-        timeValue startTime);
+        TimeValueType startTime);
 
     ~simulationFieldsDataBase() override = default;
 

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -178,4 +178,26 @@ bool PostprocessOperationAverage::write(const fileSystem &parDir) const
     return true;
 }
 
+bool PostprocessOperationAverage::write(iOstream &os) const
+{
+    if(! postprocessOperation::write(os))
+    {
+        return false;
+    }
+    if(!calculateFluctuation2_())
+    {
+        return true;
+    }
+    
+    return 
+        std::visit
+        (
+            [&](auto&& arg)->bool
+            {
+                return arg.writeFieldToVtk(os);
+            },
+            fluctuation2FieldPtr_()
+        );
+}
+
 } // namespace pFlow::postprocessData

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

@@ -195,6 +195,8 @@ public:
     /// write to os stream
     bool write(const fileSystem &parDir)const override;
 
+    bool write(iOstream& os)const override;
+
 
     /// @brief Execute average operation on field values
     /// @param weights Weight factors for particles

src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp

@@ -142,7 +142,7 @@ regionField<T> executeFluctuation2Operation
 )
 {
     const auto& regPoints = fieldAvg.regPoints();
-    regionField<T> processedField(regFieldName, regPoints, T{});
+    regionField<T> processedField(regFieldName+"_fluctuation2", regPoints, T{});
     auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)

src/PostprocessData/operation/includeMask/IncludeMask.hpp

@@ -86,12 +86,12 @@ private:
     word                    fieldName_;
     
     /// Timestamp when mask was last updated (-1 indicates never updated)
-    timeValue               lastUpdated_ = -1;
+    TimeValueType               lastUpdated_ = -1;
 
     /// Updates the mask based on current field values if needed, returns true if successful
     bool updateMask()
     {
-        timeValue t = database().currentTime();
+        TimeValueType t = database().currentTime();
         
         if( equal( t, lastUpdated_)) return true;
         
@@ -206,11 +206,11 @@ private:
 	
     std::vector<bool>       mask_;
 
-    timeValue               lastUpdated_ = -1;
+    TimeValueType               lastUpdated_ = -1;
 
     bool updateMask()
     {
-        timeValue t = database().currentTime();
+        TimeValueType t = database().currentTime();
         
         if( equal( t, lastUpdated_)) return true;
         

src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp

@@ -124,6 +124,25 @@ bool postprocessOperation::write(const fileSystem &parDir) const
     return true;
 }
 
+bool postprocessOperation::write(iOstream& os)const
+{
+    if(!regPoints().writeToSameTimeFile())
+    {
+        const auto& field = processedField();
+
+        return 
+        std::visit
+        (
+            [&](auto&& arg)->bool
+            {
+                return arg.writeFieldToVtk(os);
+            },
+            field
+        );
+    }
+    return false;
+}
+
 uniquePtr<postprocessOperation> postprocessOperation::create
 (
     const dictionary &opDict,

src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp

@@ -255,7 +255,7 @@ public:
     /// write the result to output stream (possibly a file)
     /// @param os Output stream to write the result.
     virtual 
-    bool write(iOstream& os)const {return true;}
+    bool write(iOstream& os)const;
 
     /// Create the polymorphic object using the virtual constructor.
     /// @param opDict Dictionary containing operation-specific parameters.

src/PostprocessData/operation/postprocessOperation/postprocessOperationFunctions.hpp

@@ -46,7 +46,7 @@ inline
 bool writeField
 (
     iOstream&             os,
-    timeValue             t,
+    TimeValueType             t,
     const regionField<T>  field,
     uint32                threshold,
     const T&              defValue=T{}    

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -152,11 +152,25 @@ bool pFlow::postprocessData::PostprocessComponent<RegionType, ProcessMethodType>
     }
     else
     {
-        notImplementedFunction;
+        word chNum = real2FixedStripZeros(database().time().currentTime() *1000000, 0);
+        fileSystem  file = parDir + (name() +"-"+chNum+".vtk");
+
+        auto osPtr = makeUnique<oFstream>(file);
+
+        regPoints().write(osPtr());
+
+        for(auto& operation:operatios_)
+        {
+            if(!operation->write(osPtr()))
+            {
+                fatalErrorInFunction
+                    <<"Error occurred in writing operation defined in dict "
+                    << operation->operationDict()
+                    <<endl;
                 return false;
             }
-    
-
+        }   
+    }
 
     return true;
 }