corrections for readme.md file postprocessing

postprocessPhasicFlow is now updated with new postprocessData auxFunc…

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/interaction

@@ -19,7 +19,7 @@ contactSearch
 {
     method         NBS;
 
-    updateInterval   10;
+    updateInterval   20;
 
     sizeRatio       1.1;
 

benchmarks/rotatingDrum/rotatingDrum_1mParticles/caseSetup/shapes

@@ -3,10 +3,13 @@
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
 
-objectName     particleInsertion;
+objectName     shapes;
-objectType     dicrionary;
+objectType     dictionary;
 fileFormat     ASCII;
+/*---------------------------------------------------------------------------*/
 
-active  No;          // is insertion active?
+names          (glassBead);    // names of shapes 
 
+diameters      (0.004);        // diameter of shapes 
 
+materials      (glassMat);     // material names for shapes 

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/domainDict

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

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/geometryDict

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

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/particlesDict

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

benchmarks/rotatingDrum/rotatingDrum_1mParticles/settings/settingsDict

@@ -6,13 +6,13 @@ objectName    settingsDict;
 objectType    dictionary;
 fileFormat    ASCII;
 /*---------------------------------------------------------------------------*/
-run             rotatingDrum_4MParticles;
+run             rotatingDrum_1mParticles;
 
 dt              0.00001;      // time step for integration (s)
 
 startTime       0;            // start time for simulation 
 
-endTime         10;           // end time for simulation 
+endTime         4;           // end time for simulation 
 
 saveInterval    0.2;          // time interval for saving the simulation
 
@@ -31,4 +31,4 @@ writeFormat               binary;              // data writting format (ascii or
 
 timersReport              Yes;      
 
-timersReportInterval      0.05;    
+timersReportInterval      0.01;    

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/interaction

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/caseSetup/shapes

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/cleanThisCase

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/domainDict

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/geometryDict

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/particlesDict

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

benchmarks/rotatingDrum/rotatingDrum_250kParticles/settings/settingsDict

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/caseSetup/interaction

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/cleanThisCase

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/domainDict

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/geometryDict

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/particlesDict

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

benchmarks/rotatingDrum/rotatingDrum_2mParticles/settings/settingsDict

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/interaction

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/caseSetup/shapes

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/cleanThisCase

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/domainDict

@@ -8,12 +8,16 @@ fileFormat     ASCII;
 /*---------------------------------------------------------------------------*/
 globalBox                                        // Simulation domain: every particles that goes outside this domain will be deleted
 {
-    min (-0.2 -0.2 -0.0);
+    min (-0.2 -0.2  0.0);
     max ( 0.2  0.2  1.6);
 }
 
 boundaries
 {
+    neighborListUpdateInterval 200;
+    
+    updateInterval              20;
+    
     left
     {
         type    exit;                            // other options: periodic, reflective 

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/particlesDict

@@ -29,7 +29,7 @@ positionParticles
     {
         diameter   0.003;           // minimum space between centers of particles
         
-        numPoints  4000000;         // number of particles in the simulation 
+        numPoints  4000000;          // number of particles in the simulation 
         
         axisOrder  (z x y);         // axis order for filling the space with particles
     }

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/settingsDict

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/interaction

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/caseSetup/shapes

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/cleanThisCase

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/domainDict

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/geometryDict

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/particlesDict

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

benchmarks/rotatingDrum/rotatingDrum_500kParticles/settings/settingsDict

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/interaction

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/caseSetup/shapes

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/cleanThisCase

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/runThisCase

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/domainDict

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/geometryDict

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/particlesDict

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

benchmarks/rotatingDrum/rotatingDrum_8mParticles/settings/settingsDict

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

cmake/autoComplete

@@ -1,48 +1,89 @@
 PF_cFlags="--description --help --version"
 AllTimeFolders=
 __getAllTime(){
-	files=( $(ls) )
+    # Initialize empty array for time folders
-	deleteFiles=(settings caseSetup cleanThisCase VTK runThisCase stl postprocess postProcess) 
+    local time_folders=()
-	declare -A delk
+    
-	for del in "${deleteFiles[@]}" ; do delk[$del]=1 ; done
+    # Loop through all directories in current folder
-		# Tag items to remove, based on
+    for dir in */; do
-	for k in "${!files[@]}" ; do
+        # Remove trailing slash
-		[ "${delk[${files[$k]}]-}" ] && unset 'files[k]'
+        dir=${dir%/}
-	done
+        
-		# Compaction
+        # Check if directory name is a valid floating point number
-	COMPREPLY=("${files[@]}")
+        # This pattern matches integers and floating point numbers
-	AllTimeFolders="${files[@]}"
+        if [[ $dir =~ ^[0-9]+(\.[0-9]+)?$ ]]; then
+            time_folders+=("$dir")
+        fi
+    done
+    
+    # Set completion reply to the time folders
+    COMPREPLY=("${time_folders[@]}")
+    AllTimeFolders="${time_folders[@]}"
 }
 
 __getFields(){
-	__getAllTime
+    __getAllTime
-	local -A unique_files=()
+    local -A unique_files=()
+    # Files to exclude from suggestions
+    local exclude_files=("shapeHash" "pStructure" "particleInsertion" "p" "alpha" "U" "Sp" "Su" "phi")
+    declare -A exclude_dict
+    
+    # Build exclude dictionary for faster lookups
+    for file in "${exclude_files[@]}"; do
+        exclude_dict["$file"]=1
+    done
 
-	for dir in $AllTimeFolders; do
+    for dir in $AllTimeFolders; do
-	  # Check if the directory exists
+        # Skip if not a directory
-	  if [ ! -d "$dir" ]; then
+        [ ! -d "$dir" ] && continue
-		continue  # Skip to the next directory
+        
-	  fi
+        # Find all files in this time directory
-		  
+        while IFS= read -r filename; do
-	  files_in_dir=$(find "$dir" -maxdepth 1 -type f -printf '%f\n' | sort -u)
+            # Skip empty lines and excluded files
-	  
+            [ -z "$filename" ] || [ "${exclude_dict[$filename]+exists}" ] && continue
-	# Add filenames to the associative array (duplicates will be overwritten)
+            
-	  while IFS= read -r filename; do
+            # Add to unique files
-		unique_files["$filename"]=1 # Just the key is important, value can be anything
+            unique_files["$filename"]=1
-	  done <<< "$files_in_dir"
+        done < <(find "$dir" -maxdepth 1 -type f -printf '%f\n')
-	  
+    done
-	done
+    
-	COMPREPLY=("${!unique_files[@]}")
+    # Set completion reply to the unique filenames
-	AllTimeFolders=
+    COMPREPLY=(${!unique_files[@]})
+    
+    # Clear global variable
+    AllTimeFolders=
 }
 
 _pFlowToVTK(){
-    if [ "$3" == "--time" ]  ; then
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    
+    # Check if we're completing a field
+    local is_field=0
+    for ((i=1; i<COMP_CWORD; i++)); do
+        if [[ "${COMP_WORDS[i]}" == "--fields" ]]; then
+            is_field=1
+            break
+        fi
+    done
+    
+    if [ "$prev" == "--time" ]; then
         __getAllTime
-    elif [ "$3" == "--fields" ]; then 
+    elif [ "$prev" == "--fields" ] || [ $is_field -eq 1 ]; then 
-    	__getFields
+        # We're completing field names
+        __getFields
+        # Filter the results based on the current word prefix
+        if [ -n "$cur" ]; then
+            local filtered=()
+            for item in "${COMPREPLY[@]}"; do
+                if [[ "$item" == "$cur"* ]]; then
+                    filtered+=("$item")
+                fi
+            done
+            COMPREPLY=("${filtered[@]}")
+        fi
     else
-        COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$2") )
+        COMPREPLY=( $(compgen -W "$PF_cFlags --binary --no-geometry --no-particles --out-folder --time --separate-surfaces --fields" -- "$cur") )
     fi
 }
 complete -F _pFlowToVTK pFlowToVTK

src/CMakeLists.txt

@@ -13,3 +13,5 @@ add_subdirectory(Interaction)
 
 add_subdirectory(MotionModel)
 
+add_subdirectory(PostprocessData)
+

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

@@ -97,10 +97,11 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 	const word& baseName,
 	pointStructure& pStruct,
 	const word& method,
-	const realx3Field_D& initialValField
+	const realx3Field_D& initialValField,
+	bool  keepHistory
 )
 :
-	integration(baseName, pStruct, method, initialValField),
+	integration(baseName, pStruct, method, initialValField, keepHistory),
 	realx3PointField_D
 	(
 		objectFile
@@ -108,7 +109,7 @@ pFlow::AdamsBashforth2::AdamsBashforth2
 			groupNames(baseName,"dy1"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -81,7 +81,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Destructor 
 		~AdamsBashforth2()override = default;

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -109,10 +109,11 @@ pFlow::AdamsBashforth3::AdamsBashforth3
 	const word& baseName,
 	pointStructure& pStruct,
 	const word& method,
-	const realx3Field_D& initialValField
+	const realx3Field_D& initialValField,
+	bool  keepHistory
 )
 :
-	AdamsBashforth2(baseName, pStruct, method, initialValField),
+	AdamsBashforth2(baseName, pStruct, method, initialValField, keepHistory),
 	dy2_
 	(
 		objectFile
@@ -120,7 +121,7 @@ pFlow::AdamsBashforth3::AdamsBashforth3
 			groupNames(baseName,"dy2"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory ? objectFile::WRITE_ALWAYS : objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -71,7 +71,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
 
 		/// Destructor 

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -115,10 +115,11 @@ pFlow::AdamsBashforth4::AdamsBashforth4
 	const word& baseName,
 	pointStructure& pStruct,
 	const word& method,
-	const realx3Field_D& initialValField
+	const realx3Field_D& initialValField,
+	bool keepHistory
 )
 :
-	AdamsBashforth3(baseName, pStruct, method, initialValField),
+	AdamsBashforth3(baseName, pStruct, method, initialValField, keepHistory),
 	dy3_
 	(
 		objectFile
@@ -126,7 +127,7 @@ pFlow::AdamsBashforth4::AdamsBashforth4
 			groupNames(baseName,"dy3"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -69,7 +69,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
 		
 

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -123,10 +123,11 @@ pFlow::AdamsBashforth5::AdamsBashforth5
 	const word &baseName, 
 	pointStructure &pStruct, 
 	const word &method, 
-	const realx3Field_D &initialValField
+	const realx3Field_D &initialValField,
+	bool  keepHistory
 )
 :
-	AdamsBashforth4(baseName, pStruct, method, initialValField),
+	AdamsBashforth4(baseName, pStruct, method, initialValField, keepHistory),
 	dy4_
 	(
 		objectFile
@@ -134,7 +135,7 @@ pFlow::AdamsBashforth5::AdamsBashforth5
 			groupNames(baseName,"dy4"),
 			pStruct.time().integrationFolder(),
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			keepHistory?objectFile::WRITE_ALWAYS:objectFile::WRITE_NEVER
 		),
 		pStruct,
 		zero3,

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -68,7 +68,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 		
 		
 

src/Integration/integration/integration.cpp

@@ -51,10 +51,12 @@ pFlow::integration::integration(
     const word &baseName,
     pointStructure &pStruct,
     const word &,
-    const realx3Field_D &)
+    const realx3Field_D &,
+	bool keepHistory)
     : owner_(*pStruct.owner()),
       pStruct_(pStruct),
-      baseName_(baseName)
+      baseName_(baseName),
+      keepHistory_(keepHistory)
 {}
 
 
@@ -64,12 +66,13 @@ pFlow::uniquePtr<pFlow::integration>
 		const word& baseName,
 		pointStructure& pStruct,
 		const word& method,
-		const realx3Field_D& initialValField
+		const realx3Field_D& initialValField,
+		bool  keepHistory
 )
 {
 	if( wordvCtorSelector_.search(method) )
 	{
-		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField);
+		return wordvCtorSelector_[method] (baseName, pStruct, method, initialValField, keepHistory);
 	}
 	else
 	{

src/Integration/integration/integration.hpp

@@ -24,6 +24,7 @@ Licence:
 
 #include "virtualConstructor.hpp"
 #include "pointFields.hpp"
+#include "Logical.hpp"
 
 
 namespace pFlow
@@ -63,6 +64,8 @@ private:
 		/// The base name for integration 
 		const word baseName_;
 
+		bool 	keepHistory_;
+
 protected:
 
 	bool insertValues(
@@ -83,7 +86,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 		/// Copy constructor 
 		integration(const integration&) = default;
@@ -109,9 +113,10 @@ public:
 				const word& baseName,
 				pointStructure& pStruct,
 				const word& method,
-				const realx3Field_D& initialValField
+				const realx3Field_D& initialValField,
+				bool keepHistory
 			),
-			(baseName, pStruct, method, initialValField)
+			(baseName, pStruct, method, initialValField, keepHistory)
 		);
 
 
@@ -138,6 +143,11 @@ public:
 			return owner_;
 		}
 
+		bool keepHistory()const
+		{
+			return keepHistory_;
+		}
+
 		virtual 
 		void updateBoundariesSlaveToMasterIfRequested() = 0;
 		/// return integration method 
@@ -164,7 +174,8 @@ public:
 			const word& baseName,
 			pointStructure& pStruct,
 			const word& method,
-			const realx3Field_D& initialValField);
+			const realx3Field_D& initialValField,
+			bool  keepHistory);
 
 }; // integration
 

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

@@ -248,7 +248,8 @@ pFlow::grainParticles::grainParticles(
 			"rVelocity",
 			dynPointStruct(),
 			intMethod,
-			rAcceleration_.field()
+			rAcceleration_.field(),
+			control.keepIntegrationHistory()
 		);
 
 	if( !rVelIntegration_ )

src/Particles/GrainParticles/grainShape/grainShape.cpp

@@ -73,6 +73,18 @@ pFlow::grainShape::grainShape
 	}
 }
 
+pFlow::grainShape::grainShape
+(
+	const word &shapeType, 
+	const word &fileName, 
+	repository *owner, 
+	const property &prop
+)
+:
+	grainShape(fileName, owner, prop)
+{
+}
+
 pFlow::real pFlow::grainShape::maxBoundingSphere() const
 {
     return max(grainDiameters_);
@@ -99,9 +111,12 @@ pFlow::real pFlow::grainShape::boundingDiameter(uint32 index) const
 	{
 		return grainDiameters_[index];
 	}
-	fatalErrorInFunction<<"Invalid index for diameter "<<
+
-	index<<endl;
+	fatalErrorInFunction
+		<<"Invalid index for diameter "
+		<<index<<endl;
 	fatalExit;
+	
 	return 0.0;
 }
 
@@ -122,13 +137,17 @@ pFlow::real pFlow::grainShape::coarseGrainFactor(uint32 index) const
 	return 0.0;
 }
 
+pFlow::realVector pFlow::grainShape::volume() const
+{
+    return realVector("volume", Pi/6*pow(grainDiameters_,(real)3.0));
+}
 
 pFlow::realVector pFlow::grainShape::coarseGrainFactor() const
 {
     return coarseGrainFactor_;
 }
 
-pFlow::real pFlow::grainShape::orginalDiameter(uint32 index) const
+pFlow::real pFlow::grainShape::originalDiameter(uint32 index) const
 {
 	if(indexValid(index))
 	{
@@ -142,7 +161,7 @@ pFlow::real pFlow::grainShape::orginalDiameter(uint32 index) const
 
 
 
-pFlow::realVector pFlow::grainShape::orginalDiameter() const
+pFlow::realVector pFlow::grainShape::originalDiameter() const
 {
     return sphereDiameters_;
 }

src/Particles/GrainParticles/grainShape/grainShape.hpp

@@ -32,9 +32,13 @@ class grainShape
 {
 private:
 	
-	// - diameter of spheres
+	/// diameter of grains
 	realVector 			grainDiameters_;
+
+	/// diameter of spheres
 	realVector 			sphereDiameters_;
+	
+	/// course-grain factor 
 	realVector 			coarseGrainFactor_;
 
 	
@@ -54,9 +58,21 @@ public:
 		repository* owner,
 		const property& prop);
 
+	grainShape(
+		const word& shapeType,
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
 	
 	~grainShape() override = default;
 
+	add_vCtor
+	(
+		shape, 
+		grainShape,
+		word
+	);
 	//// - Methods
 
 	real maxBoundingSphere()const override;
@@ -68,14 +84,16 @@ public:
 	real boundingDiameter(uint32 index)const override;
 
 	realVector boundingDiameter()const override;
+	
+	realVector volume()const override;
 
-		real coarseGrainFactor(uint32 index)const ;
+	real coarseGrainFactor(uint32 index)const ;
 
 	realVector coarseGrainFactor()const ;
 
-	real orginalDiameter(uint32 index)const ;
+	real originalDiameter(uint32 index)const ;
 
-	realVector orginalDiameter()const ;
+	realVector originalDiameter()const ;
 
 	bool mass(uint32 index, real& m)const override;
 

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -229,7 +229,8 @@ pFlow::sphereParticles::sphereParticles(
 			"rVelocity",
 			dynPointStruct(),
 			intMethod,
-			rAcceleration_.field()
+			rAcceleration_.field(),
+			control.keepIntegrationHistory()
 		);
 
 	if( !rVelIntegration_ )

src/Particles/SphereParticles/sphereShape/sphereShape.cpp

@@ -68,6 +68,18 @@ pFlow::sphereShape::sphereShape
 	}
 }
 
+pFlow::sphereShape::sphereShape
+(
+	const word &shapeType, 
+	const word &fileName, 
+	repository *owner, 
+	const property &prop
+)
+:
+	sphereShape(fileName, owner, prop)
+{
+}
+
 pFlow::real pFlow::sphereShape::maxBoundingSphere() const
 {
     return max(diameters_);
@@ -105,6 +117,11 @@ pFlow::realVector pFlow::sphereShape::boundingDiameter() const
     return diameters_;
 }
 
+pFlow::realVector pFlow::sphereShape::volume() const
+{
+    return realVector("volume", Pi/6*pow(diameters_,(real)3.0));
+}
+
 bool pFlow::sphereShape::mass(uint32 index, real &m) const
 {
     if( indexValid(index) )

src/Particles/SphereParticles/sphereShape/sphereShape.hpp

@@ -51,9 +51,22 @@ public:
 		repository* owner,
 		const property& prop);
 
+	sphereShape(
+		const word& shapeType,
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
 	
 	~sphereShape() override = default;
 
+	add_vCtor
+	(
+		shape,
+		sphereShape,
+		word
+	);
+
 	//// - Methods
 
 	real maxBoundingSphere()const override;
@@ -66,6 +79,8 @@ public:
 
 	realVector boundingDiameter()const override;
 
+	realVector volume()const override;
+
 	bool mass(uint32 index, real& m)const override;
 
 	real mass(uint32 index) const override;

src/Particles/dynamicPointStructure/dynamicPointStructure.cpp

@@ -58,13 +58,14 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 {
 	REPORT(1)<< "Creating integration method "<<
 		Green_Text(integrationMethod_)<<" for dynamicPointStructure."<<END_REPORT;
-
+	
 	integrationPos_ = integration::create
 	(
 		"pStructPosition",
 		*this,
 		integrationMethod_,
-		velocity_.field()
+		velocity_.field(),
+		control.keepIntegrationHistory()
 	);
 
 	if( !integrationPos_ )
@@ -79,7 +80,8 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		"pStructVelocity",
 		*this,
 		integrationMethod_,
-		acceleration_.field()
+		acceleration_.field(),
+		control.keepIntegrationHistory()
 	);
 
 	if( !integrationVel_ )

src/Particles/particles/particles.hpp

@@ -186,7 +186,7 @@ public:
 	}
 
 	inline 
-	uint maxId()const
+	uint32 maxId()const
 	{
 		return idHandler_().maxId();
 	}

src/Particles/particles/shape/shape.cpp

@@ -1,5 +1,24 @@
-#include "shape.hpp"
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
+-----------------------------------------------------------------------------*/
+
+#include "shape.hpp"
 
 bool pFlow::shape::findPropertyIds()
 {	
@@ -62,6 +81,18 @@ pFlow::shape::shape
 
 }
 
+pFlow::shape::shape
+(
+	const word &shapeType, 
+	const word &fileName, 
+	repository *owner, 
+	const property &prop
+)
+:
+	shape(fileName, owner, prop)
+{
+}
+
 bool pFlow::shape::writeToDict(dictionary &dict)const
 {
 	if(!baseShapeNames::writeToDict(dict))return false;
@@ -74,4 +105,37 @@ bool pFlow::shape::writeToDict(dictionary &dict)const
 	}
 
 	return true;
-}
+}
+
+pFlow::uniquePtr<pFlow::shape> pFlow::shape::create
+(
+	const word &shapeType, 
+	const word &fileName, 
+	repository *owner, 
+	const property &prop
+)
+{
+	word type = angleBracketsNames("shape", shapeType);
+
+	if( wordvCtorSelector_.search(type) )
+	{
+		auto objPtr = 
+		wordvCtorSelector_[type]
+			(shapeType, fileName, owner, prop);
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< 
+			type << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		fatalExit;
+		return nullptr;
+	}
+
+}

src/Particles/particles/shape/shape.hpp

@@ -60,9 +60,28 @@ public:
 		const word& fileName,
 		repository* owner,
 		const property& prop);
+	
+	shape(
+		const word& shapeType,
+		const word& fileName,
+		repository* owner,
+		const property& prop);
 
     ~shape() override=default;
 
+	create_vCtor
+	(
+		shape,
+		word,
+		(
+			const word& shapeType,
+			const word& fileName,
+			repository* owner,
+			const property& prop
+		),
+		(shapeType, fileName, owner, prop)	
+	);
+
 	inline 
 	const auto& properties()const
 	{
@@ -148,6 +167,9 @@ public:
 	virtual
 	realVector boundingDiameter()const = 0;
 
+	virtual 
+	realVector volume()const = 0;
+
     virtual 
     bool mass(uint32 index, real& m)const = 0;
 
@@ -187,6 +209,13 @@ public:
     virtual 
     real Inertial_zz(uint32 index)const = 0;
 
+	static
+	uniquePtr<shape> create(
+		const word& shapeType,
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
 };
 
 }

src/PostprocessData/CMakeLists.txt

@@ -4,12 +4,14 @@ set(SourceFiles
 
     # Fields database
     fieldsDataBase/fieldsDataBase.cpp
+    fieldsDataBase/simulationFieldsDataBase.cpp
   
     # Regions
     region/regionPoints/regionPoints/regionPoints.cpp
     region/regionPoints/sphereRegionPoints/sphereRegionPoints.cpp
     region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
     region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
+    region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
 
     # Postprocess components
     postprocessComponent/postprocessComponent/postprocessComponent.cpp
@@ -19,6 +21,9 @@ set(SourceFiles
     # Operations
     operation/postprocessOperation/postprocessOperation.cpp
     operation/PostprocessOperation/PostprocessOperationSum.cpp
+    operation/PostprocessOperation/PostprocessOperationAverage.cpp
+    operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp
+    
     operation/includeMask/includeMask.cpp
     operation/includeMask/IncludeMasks.cpp
 

src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp

@@ -21,25 +21,38 @@ Licence:
 #include <regex>
 
 #include "vocabs.hpp"
-#include "Time.hpp"
+#include "systemControl.hpp"
 #include "fieldsDataBase.hpp"
 #include "fieldFunctions.hpp"
-#include "dynamicPointStructure.hpp"
+#include "dictionary.hpp"
 
-bool pFlow::fieldsDataBase::checkForUpdate(const word &name, bool forceUpdate)
+namespace pFlow
+{
+
+bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+
+}
+
+bool pFlow::fieldsDataBase::loadPointStructureToTime()
+{
+    return false;
+}
+
+
+bool pFlow::fieldsDataBase::checkForUpdate(const word &compoundName, bool forceUpdate)
 {
     auto t = currentTime();
     bool shouldUpdate = false;
     
-    if(auto [iter, found]= captureTime_.findIf(name); found)
+    if(auto [iter, found]= captureTime_.findIf(compoundName); found)
     {
-        shouldUpdate = iter->second < t;
+        shouldUpdate = iter->second < t || forceUpdate;
         iter->second = t;
     }
     else
     {
         shouldUpdate = true;
-        captureTime_.insertIf(name, t);
+        captureTime_.insertIf(compoundName, t);
     }
 
     return shouldUpdate;
@@ -70,6 +83,177 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetRealField(const word
         field.size());
 }
 
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetVolume(bool forceUpdate)
+{
+    const word fName = "volume";
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+        
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        auto vols = getShape().volume();
+        
+        FieldTypeHost<real> volField
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< volField.size(); i++)
+        {
+            volField[i] = vols[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(volField)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetDensity(bool forceUpdate)
+{
+    const word fName = "density";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        const auto dens = getShape().density();
+
+        FieldTypeHost<real> denFeild
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< denFeild.size(); i++)
+        {
+            denFeild[i] = dens[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(denFeild)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetOne(bool forceUpdate)
+{
+    const word fName = "one";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            FieldTypeHost<real> 
+            (
+                fName, 
+                "value",
+                pointFieldSize(),
+                1.0
+            )  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetMass(bool forceUpdate)
+{
+    const word fName = "mass";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        const auto ms = getShape().mass();
+
+        FieldTypeHost<real> massField
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< massField.size(); i++)
+        {
+            massField[i] = ms[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(massField)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetI(bool forceUpdate)
+{
+    const word fName = "I";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        const auto Is = getShape().Inertia();
+
+        FieldTypeHost<real> IField
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< IField.size(); i++)
+        {
+            IField[i] = Is[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(IField)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
 bool pFlow::fieldsDataBase::findFunction(
     const word &compoundFieldName,
     word &fieldName,
@@ -274,95 +458,145 @@ bool pFlow::fieldsDataBase::inputOutputType
     return false;
 }
 
-pFlow::fieldsDataBase::fieldsDataBase(const Time &time, bool inSimulation)
+pFlow::fieldsDataBase::fieldsDataBase
-: 
+(
-    time_(time), 
+    systemControl& control,
+    const dictionary& postDict,
+    bool inSimulation,
+    timeValue startTime
+)
+:
+    time_(control.time()),
     inSimulation_(inSimulation)
-{}
+{
+    if(!inSimulation_)
+    {
+        shapeType_ = postDict.getValOrSet<word>("shapeType", "");
+
+        if(shapeType_.empty())
+        {
+            WARNING
+            << "shapeType is not set in dictionary: "
+            << postDict.globalName()
+            << " and you may need to set for some postprocess operations"<<END_WARNING;
+        }
+    }
+    else
+    {
+        // this is not required for execution during simulation 
+        shapeType_ = "";
+    }
+}
 
 pFlow::timeValue pFlow::fieldsDataBase::currentTime() const
 {
     return time_.currentTime();
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldTypeNameFunction
 (
-    const word& compoundName, 
+    const word& compoundName,
-    word& fieldName,
+    word& pointFieldName, 
     word& originalType, 
     word& typeAfterFunction,
     Functions& func
-)
+)const
 {
-    if( !findFunction(compoundName, fieldName, func))
+    if( !findFunction(compoundName, pointFieldName, func))
     {
-        fatalErrorInFunction;
+        fatalErrorInFunction
+            <<"Error in processing function for field: "
+            <<compoundName<<endl;;
         return false;
     }
 
-    if(reservedFieldNames_.contains(fieldName))
+    if( reservedFieldNames_.contains(pointFieldName))
     {
-        originalType = reservedFieldNames_.find(fieldName)->second;
+        // The name is in the reserved fields 
+        originalType = reservedFieldNames_.find(pointFieldName)->second;
     }
     else
     {
-        word fieldTypeName =  time_.lookupObjectTypeName(fieldName);
+        // the name is in the Time object 
-        word space;
+        if(pointFieldNameExists(pointFieldName))
-        if(!pointFieldGetType(fieldTypeName, originalType, space))
         {
-            fatalErrorInFunction<<"Cannot extract type name from "<< fieldTypeName<<endl;
+            originalType = getPointFieldType(pointFieldName);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "The field name: "<< pointFieldName
+                << " is not in the Time object.\n"
+                <<" Avaiable ones are: \n"<< time().objectNames()<<endl;
             return false;
         }
     }
-    
+
-    word outputType;
+    if(!inputOutputType(func, originalType, typeAfterFunction))
-    if(!inputOutputType(func, originalType, outputType))
     {
-        fatalErrorInFunction<<"Cannnot determine the input and output types for: "<< compoundName<<endl;
+        fatalErrorInFunction
+            <<"Cannnot determine the input and output types for: "
+            << compoundName<<endl;
         return false;
     }
-    typeAfterFunction = outputType;
+
     return true;
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldType
 (
-    const word &compoundName, 
+    const word& compoundName,
-    word &originalType, 
+    word& originalType, 
-    word &typeAfterFunction
+    word& typeAfterFunction
-)
+) const
 {
     Functions func;
     word fieldName;
-    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldType
 (
     const word &compoundName, 
     word &typeAfterFunction
-)
+) const
 {
     Functions func;
     word originalType;
     word fieldName;
-    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
 }
 
 pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updatePoints(bool forceUpdate)
 {
-    
+    const word fName = "position";
-    bool shouldUpdate = checkForUpdate("position", forceUpdate);
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
         
     if(shouldUpdate)
     {
+        // load pointStructure 
+        if(!loadPointStructureToTime())
+        {
+            fatalErrorInFunction<< "Error in loading pointStructure to Time object."<<endl;
+            fatalExit;
+        }
         const auto& pstruct = pStruct(); 
-        allFields_.emplaceBackOrReplace<PointsTypeHost>(
+        allFields_.emplaceBackOrReplace<PointsTypeHost>
-            "position", 
+        (
-            pstruct.activePointsHost());
+            fName, 
+            pstruct.activePointsHost()
+        );
     }
    
-    auto& points = allFields_.getObject<PointsTypeHost>("position");
+    auto& points = allFields_.getObject<PointsTypeHost>(fName);
     
     return span<realx3>(
         points.data(), 
@@ -379,7 +613,12 @@ pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updateFieldRealx3
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func) )
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -410,7 +649,12 @@ pFlow::span<pFlow::realx4> pFlow::fieldsDataBase::updateFieldRealx4
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func))
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -441,7 +685,12 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::updateFieldReal
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func) )
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -609,11 +858,9 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
     bool forceUpdate
 )
 {
-
+    word typeAfterFunction, originalType;
-    word originalType, typeAfterFunction, fieldName;
-    Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldType(compoundName, originalType, typeAfterFunction))
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -621,8 +868,7 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
         return span<real>(nullptr, 0);
     }
 
-
+    if(  typeAfterFunction == getTypeName<realx3>()  )
-    if(  typeAfterFunction== getTypeName<realx3>()  )
     {
         return updateFieldRealx3(compoundName, forceUpdate);
     }
@@ -642,19 +888,47 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
     }
 }
 
-const pFlow::pointStructure &pFlow::fieldsDataBase::pStruct() const
+
+
+pFlow::uniquePtr<pFlow::fieldsDataBase> 
+    pFlow::fieldsDataBase::create
+(
+    systemControl& control, 
+    const dictionary& postDict, 
+    bool inSimulation,
+    timeValue startTime
+)
 {
-    if(inSimulation_)
+    word dbType;
+    if(inSimulation)
     {
-        return 
+        dbType = "fieldsDataBase<simulation>";
-            static_cast<const pointStructure&>(
-                time_.lookupObject<dynamicPointStructure>(pointStructureFile__)
-            );
     }
     else
     {
-        return time_.lookupObject<pointStructure>(pointStructureFile__);
+        dbType = "fieldsDataBase<postSimulation>";
     }
+    
+    if( boolvCtorSelector_.search(dbType) )
+    {
+        auto objPtr = 
+            boolvCtorSelector_[dbType](control, postDict, inSimulation, startTime);
+        return objPtr;
+    }
+    else
+    {
+        printKeys
+        ( 
+            fatalError << "Ctor Selector "<< 
+            dbType << " does not exist. \n"
+            <<"Available ones are: \n\n"
+            ,
+            boolvCtorSelector_
+        );
+        fatalExit;
+    }
+    
+    return nullptr;
 }
 
 bool pFlow::pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace)

src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp

@@ -21,52 +21,24 @@ Licence:
 #ifndef __fieldsDataBase_hpp__
 #define __fieldsDataBase_hpp__
 
-#include <variant>
+#include "fieldsDataBaseDclr.hpp"
-#include <concepts>
-
 
 #include "pointStructure.hpp"
 #include "pointFields.hpp"
 #include "anyList.hpp"
 #include "Map.hpp"
-#include "span.hpp"
+#include "shape.hpp"
+
 
 namespace pFlow
 {
 
+class dictionary;
+class systemControl;
 class Time;
 
-bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
-
-template<typename T>
-concept ValidFieldType = 
-    std::same_as<T, real> || 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4> || 
-    std::same_as<T, uint32>;
-
-template<typename T> 
-concept VectorType = 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4>;
-
-template<typename T>
-concept ScalarType = 
-    std::same_as<T, real>;
-
-
-template<typename T>
-concept ValidRegionFieldType = 
-    std::same_as<T, real> || 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4> ;
-
-
-using allPointFieldTypes = std::variant<span<real>, span<realx3>, span<realx4>>; 
-
 
 class fieldsDataBase
-
 {
 private:
 
@@ -99,9 +71,6 @@ private:
 
 
     // - Member variables
-    
-        /// const reference to the Time object
-        const Time&         time_;
 
         /// List to store all the point fields 
         anyList             allFields_;
@@ -109,114 +78,105 @@ private:
         /// Map to store the last capture time of each field
         wordMap<timeValue>  captureTime_;
 
-        /// 
+        /// Reference to the Time object
-        bool                inSimulation_ = false;
+        Time&               time_; 
-        
-        static const inline std::map<word, word> reservedFieldNames_
-        {
-            {"position", "realx3"},
-            {"one", "real"}
-        };
-        
-    // - Private member functions
 
+        /// Flag indicating if we're in simulation mode
+        bool                inSimulation_ = false;
+
+        word                shapeType_;
+
+protected:
+
+        /// check if pointField name exists in Time or time folder 
+        virtual 
+        bool pointFieldNameExists(const word& name)const = 0;
+        
+        /// Loads a pointField with a given name to the Time object 
+        virtual 
+        bool loadPointFieldToTime(const word& name)= 0;
+
+        virtual 
+        bool loadPointStructureToTime()=0;
+        
+        virtual 
+        const shape& getShape() const= 0;
+
+        const word& shapeTypeName()const
+        {
+            return shapeType_;
+        }
+
+        /// get the type name of the pointField in the Time object 
+        virtual 
+        word getPointFieldType(const word& name)const = 0;
+
+        /// Checks if a field needs to be updated based on capture time
+        bool checkForUpdate(const word& compoundName, bool forceUpdate = false);
+    
+        /// @brief return the size of pointStructure 
         uint32 pointFieldSize()
         {
             auto s = updatePoints();
             return s.size();
         }
-        
-        /// @brief Checks if a field needs to be updated.
-        /// @param name The name of the field to check.
-        /// @param forceUpdate Forces an update if true. Defaults to `false`.
-        /// @return `true` if the field needs updating or `forceUpdate` is true, `false` otherwise.
-        bool checkForUpdate(const word& name, bool forceUpdate = false);
 
-        /**
-         * @brief Updates and retrieves a point field of a specified type from the database.
-         * 
-         * This is a template function that updates and retrieves a point field of type `T`
-         * from the database. It checks if the field needs to be updated based on the last
-         * capture time or if a forced update is requested. If an update is necessary, it
-         * retrieves the latest data for the field.
-         * 
-         * @tparam T The type of the point field to update and retrieve. Must be a ValidFieldType.
-         * @param name The name of the field to update.
-         * @param forceUpdate A boolean flag indicating whether to force an update of the field
-         *                    regardless of its current state. Defaults to `false`.
-         * 
-         * @return A span of `T` representing the updated field data.
-         * 
-         * @throws message If the field cannot be found in the database or if there is a type mismatch.
-         */
         template<ValidFieldType T>
         span<T> updateField(const word& name, bool forceUpdate = false);
+
+        template<ValidFieldType T>
+        span<T> updateReservedField(const word& name, bool forceUpdate = false);
         
-        /// @brief Creates a new real field or retrieves an existing one.
-        /// 
-        /// If a field with the given name already exists, it returns a span to that field.
-        /// If the field does not exist, it creates a new real field with the given name
-        /// and returns a span to the newly created field.
-        /// 
-        /// @param name The name of the field to create or retrieve.
-        /// @return span<real> of the field 
         span<real> createOrGetRealField(const word& name);
 
-        /**
+        span<real> createOrGetVolume(bool forceUpdate=false);
-         * @brief Parses a compound field name to extract the base field name and function.
+
-         * 
+        span<real> createOrGetDensity(bool forceUpdate=false);
-         * This function takes a compound field name, which may include a function applied
+
-         * to a base field (e.g., "mag(velocity)"), and parses it to extract the base field
+        span<real> createOrGetOne(bool forceUpdate=false);
-         * name (e.g., "velocity") and the function to be applied (e.g., `Functions::Magnitude`).
+
-         * 
+        span<real> createOrGetMass(bool forceUpdate=false);
-         * The function supports the following syntax for compound field names:
+
-         * - `fieldName` (no function applied)
+        span<real> createOrGetI(bool forceUpdate=false);
-         * - `functionName(fieldName)`
+
-         * 
+        /// Map of reserved field names with their corresponding types
-         * Supported function names are defined in the `Functions` enum.
+        static const inline std::map<word, word> reservedFieldNames_
-         * 
+        {
-         * @param compoundFieldName The compound field name to parse.
+            {"position", "realx3"},
-         * @param fieldName A reference to a `word` where the extracted base field name
+            {"one", "real"},
-         *                  will be stored.
+            {"volume", "real"},
-         * @param func A reference to a `Functions` enum where the identified function
+            {"density", "real"},
-         *             will be stored. If no function is applied, this will be set to
+            {"mass", "real"},
-         *             `Functions::None`.
+            {"I", "real"}
-         * 
+        };
-         * @return `true` if the compound field name was successfully parsed and the base
+        
-         *         field name and function were extracted, `false` otherwise.
-         * 
-         * @note The function modifies the `fieldName` and `func` parameters to return the
-         *       extracted information.
-         */
         static 
         bool findFunction(
             const word& compoundFieldName, 
             word& fieldName, 
             fieldsDataBase::Functions& func );
         
-        /**
-         * @brief Determines the input and output types for a given function.
-         * 
-         * This function takes a `Functions` enum value and an input type as a string
-         * and determines the corresponding output type based on the function being applied.
-         * 
-         * @param func The function for which to determine the input and output types.
-         * @param inputType The input type as a string.
-         * @param outputType A reference to a string where the determined output type will be stored.
-         * 
-         * @return `true` if the input and output types were successfully determined, `false` otherwise.
-         */
         static
         bool inputOutputType(
             fieldsDataBase::Functions func, 
             const word& inputType, 
             word& outputType);
 
+    
+
 public:
     
+    // - Type info
+            
+        TypeInfo("fieldsDataBase");
+
     // - constructors 
 
-        fieldsDataBase(const Time& time, bool inSimulation);
+        fieldsDataBase(
+            systemControl& control, 
+            const dictionary& postDict,
+            bool inSimulation,
+            timeValue startTime);
 
         /// no copy constructor
         fieldsDataBase(const fieldsDataBase&) = delete;
@@ -231,140 +191,75 @@ public:
         fieldsDataBase& operator=(fieldsDataBase&&) = delete;
         
         /// destructor
-        ~fieldsDataBase() = default;
+        virtual ~fieldsDataBase() = default;
+
+        create_vCtor
+        (
+            fieldsDataBase,
+            bool,
+            (
+                systemControl& control, 
+                const dictionary& postDict, 
+                bool inSimulation,
+                timeValue startTime
+            ),
+            (control, postDict, inSimulation, startTime)
+        );
+  
     
     // - Public Access Functions
         /// returns the current time
         timeValue currentTime()const;
 
-        /// const ref to object Time
+        /// const ref to object Time 
-        const Time& time()const
+        const Time& time()const 
         {
             return time_;
         }
-    
+        
+        Time& time()
+        {
+            return time_;
+        }
+
     // - Public Member Functions
 
-        /**
+        bool getFieldTypeNameFunction
-         * @brief Retrieves the type of a point field based on its compound name.
+        (
-         * 
-         * This function attempts to extract the type of a point field from its compound name.
-         * The compound name may include additional information such as a function or operation
-         * applied to the field, ie. mag(velcoty). If the type is successfully determined, it  
-         * is stored in the provided `typeName` parameter.
-         * 
-         * @param compoundName The compound name of the field, which may include additional
-         *                     information about operations or functions applied to the field.
-         * @param originalType A reference to a `word` where the original type name is obtained.
-         *                     This will be set to the type of the field before any function is applied.
-         * @param typeAfterFunction A reference to a `word` where the type name after applying                   
-         *                          the function is obtained.
-         * @param func the applied function to the field.
-         *  
-         * @return `true` if the type was successfully determined and stored in `typeName`,
-         *         `false` otherwise.
-         */
-        bool getPointFieldType(
             const word& compoundName,
-            word& fieldName, 
+            word& pointFieldName, 
             word& originalType, 
             word& typeAfterFunction,
-            Functions& func);   
+            Functions& func
-        
+        )const;
-        /// overload for the function getPointFieldType without `func` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& originalType, 
-            word& typeAfterFunction);
-        
-        /// overload for function getPointFieldType without `originalType` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& typeAfterFunction);
 
-        /**
+        
-         * @brief Updates the points data and returns a span to the updated points.
+        bool getFieldType
-         * 
+        (
-         * This function ensures that the points data is up-to-date by checking if an update
+            const word& compoundName,
-         * is necessary. If the data is outdated or if a forced update is requested, it retrieves
+            word& originalType, 
-         * the latest points data and stores it in the internal fields database. The function
+            word& typeAfterFunction
-         * then returns a span to the updated points data for further use.
+        ) const;
-         * 
+
-         * @param forceUpdate A boolean flag indicating whether to force an update of the points
+        bool getFieldType
-         *                    data regardless of its current state. Defaults to `false`.
+        (
-         * 
+            const word& compoundName,
-         * @return A span of `realx3` representing the updated points data.
+            word& typeAfterFunction
-         */
+        ) const;
+
+        
+        /// update pointStructure if necessary 
         span<realx3> updatePoints(bool forceUpdate = false);
         
-        /**
+        /// update a field with realx3 type 
-         * @brief Updates and retrieves a realx3 point field from the database.
-         * 
-         * This function retrieves or updates a realx3 field based on its compound name.
-         * The compound name cannot include any function operation.
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated realx3 field data
-         * 
-         * @throws message If the field type is not compatible with realx3 or if the field
-         *                 cannot be found in the database
-         */
         span<realx3> updateFieldRealx3(
             const word& compoundName, 
             bool forceUpdate = false);
         
-        /**
-         * @brief Updates and retrieves a realx4 point field from the database.
-         * 
-         * This function retrieves or updates a realx4 field based on its compound name.
-         * The compound name cannot include any function operation.
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated realx3 field data
-         * 
-         * @throws message If the field type is not compatible with realx4 or if the field
-         *                 cannot be found in the database
-         */
         span<realx4> updateFieldRealx4(
             const word& compoundName, 
             bool forceUpdate = false);
         
-        /**
-         * @brief Updates and retrieves a real point field from the database.
-         * 
-         * This function retrieves or updates a real field based on its compound name.
-         * The compound name may include a function operation (e.g., abs, square, etc.).
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database and apply the specified function.
-         * 
-         * Supported functions include:
-         * - None: Retrieves the field as is.
-         * - abs: Computes the absolute value of the field.
-         * - square: Computes the square of the field.
-         * - cube: Computes the cube of the field.
-         * - sqrt: Computes the square root of the field.
-         * - mag, magSquare, magCube, magSqrt: Compute magnitude-related operations for vector fields.
-         * - component(x, y, z): Extracts a specific component from a vector field.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation.
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated real field data.
-         * 
-         * @throws message If the field type is not compatible with real or if the field
-         *                 cannot be found in the database.
-         */
         span<real> updateFieldReal(
             const word& compoundName, 
             bool forceUpdate = false);
@@ -378,11 +273,45 @@ public:
         allPointFieldTypes updateFieldAll(
             const word& compoundName, 
             bool forceUpdate = false);
+        
+        virtual 
+        const pointStructure& pStruct()const = 0;
+        
+        /// Get the next avaiable time folder after the current time folder 
+        /// This is only used for post-simulation processing 
+        virtual
+        timeValue getNextTimeFolder()const
+        {
+            return -1.0;
+        }
 
-        const pointStructure& pStruct()const;
+        /// Sets the current folder to the next time folder.
+        /// This is used only for post-simulation processing
+        /// @returns the time value of the next folder.
+        virtual 
+        timeValue setToNextTimeFolder()
+        {
+            return -1.0;
+        }
+        
+        /// Skips the next time folder.
+        /// This is used only for post-simulation processing 
+        /// @returns the time value of the skipped folder 
+        virtual 
+        timeValue skipNextTimeFolder()
+        {
+            return -1.0;
+        }
+
+    static
+    uniquePtr<fieldsDataBase> create(
+        systemControl& control, 
+        const dictionary& postDict, 
+        bool inSimulation,
+        timeValue startTime);
 };
 
-} // nampespace pFlow
+} // namespace pFlow
 
 #include "fieldsDataBaseTemplates.cpp"
 

src/PostprocessData/fieldsDataBase/fieldsDataBaseDclr.hpp

@@ -0,0 +1,45 @@
+
+#ifndef __fieldsDataBaseDclr_hpp__
+#define __fieldsDataBaseDclr_hpp__
+
+#include <variant>
+#include <concepts>
+#include <type_traits>
+
+#include "types.hpp"
+#include "span.hpp"
+
+namespace pFlow
+{
+
+
+template<typename T>
+concept ValidFieldType = 
+    std::same_as<T, real> || 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4> || 
+    std::same_as<T, uint32>;
+
+template<typename T> 
+concept VectorType = 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4>;
+
+template<typename T>
+concept ScalarType = 
+    std::same_as<T, real>;
+
+
+template<typename T>
+concept ValidRegionFieldType = 
+    std::same_as<T, real> || 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4> ;
+
+
+using allPointFieldTypes = std::variant<span<real>, span<realx3>, span<realx4>>; 
+
+
+} // namespace pFlow
+
+#endif //__fieldsDataBaseDclr_hpp__

src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp

@@ -32,40 +32,25 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
         
     if(shouldUpdate)
     {
-        if(name == "one")
+        if(reservedFieldNames_.contains(name))
-        {   
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>
-            (
-                "one",
-                FieldTypeHost<T>
-                (
-                    "one", 
-                    "value",
-                    pointFieldSize(),
-                    T{1}
-                )
-            );
-        }
-        else if( name == "position")
         {
-            if constexpr( std::same_as<T, realx3>)
+            return updateReservedField<T>(name, true);
+        }   
+        else
+        {
+            if( loadPointFieldToTime(name) )
             {
-                return updatePoints(forceUpdate);            
+                const auto& pField = time_.lookupObject<pointField_D<T>>(name);
+                allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
+                name, 
+                pField.activeValuesHost());
             }
             else
             {
-                fatalErrorInFunction<< "Error in getting the type name of field: "<<
+                fatalErrorInFunction<<"Error in loading the pointField "<<name<<endl;
-                    name<<", with type name: "<< getTypeName<T>() <<endl;
                 fatalExit;
-                return span<T>(nullptr, 0);
             }
-        }
+            
-        else
-        {
-            const auto& pField = time_.lookupObject<pointField_D<T>>(name);
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
-            name, 
-            pField.activeValuesHost());
         }
     }
 
@@ -77,4 +62,110 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
     
 }
 
+
+template<pFlow::ValidFieldType T>
+inline
+pFlow::span<T> pFlow::fieldsDataBase::updateReservedField
+(
+    const word& name, 
+    bool forceUpdate
+)
+{
+    if(name == "one")
+    {   
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetOne(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field one."<<endl;
+            fatalExit;
+        }
+    }
+    else if(name == "volume")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetVolume(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field volume."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "density")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetDensity(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field density."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "mass")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetMass(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field mass."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "I")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetI(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field I."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "position")
+    {
+        if constexpr( std::same_as<T, realx3>)
+        {
+            return updatePoints(forceUpdate);            
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field position."<<endl;
+            fatalExit;
+        }
+    }
+    
+    fatalErrorInFunction<<"Not supported field name "<<endl;
+    fatalExit;
+    return span<T>(nullptr, 0);
+
+}
+
 #endif //__fieldsDataBaseTemplates_hpp__

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.cpp

@@ -0,0 +1,71 @@
+#include "Time.hpp"
+#include "simulationFieldsDataBase.hpp"
+#include "dynamicPointStructure.hpp"
+#include "vocabs.hpp"
+
+namespace pFlow
+{
+    bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+}
+
+bool pFlow::simulationFieldsDataBase::pointFieldNameExists(const word &name) const
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::simulationFieldsDataBase::loadPointFieldToTime(const word &name)
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::simulationFieldsDataBase::loadPointStructureToTime()
+{
+    // it is already in the Time object 
+    return time().lookupObjectName(pointStructureFile__);
+}
+
+
+const pFlow::shape& pFlow::simulationFieldsDataBase::getShape() const
+{
+    return shape_;
+}
+
+pFlow::word pFlow::simulationFieldsDataBase::getPointFieldType(const word &name) const
+{
+    word pfType =  time().lookupObjectTypeName(name);
+    word type, space;
+    if(!pointFieldGetType(pfType, type, space))
+    {
+        fatalErrorInFunction
+            <<"Error in retriving the type of pointField "
+            << pfType<<endl;
+        fatalExit;
+    }
+    return type;
+}
+
+pFlow::simulationFieldsDataBase::simulationFieldsDataBase
+(
+    systemControl &control,
+    const dictionary& postDict, 
+    bool inSimulation,
+    timeValue startTime
+)
+:
+    fieldsDataBase(control, postDict, inSimulation, startTime),
+    shape_
+    (
+        dynamic_cast<const shape&>(*control.caseSetup().lookupObjectPtr(shapeFile__))
+    )
+{
+}
+
+const pFlow::pointStructure &pFlow::simulationFieldsDataBase::pStruct() const
+{
+    return 
+    static_cast<const pointStructure&>
+    (
+        time().lookupObject<dynamicPointStructure>(pointStructureFile__)
+    );
+}
+

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.hpp

@@ -0,0 +1,82 @@
+/*------------------------------- 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 __simulationFieldsDataBase_hpp__
+#define __simulationFieldsDataBase_hpp__
+
+#include "fieldsDataBase.hpp"
+
+namespace pFlow
+{
+
+class simulationFieldsDataBase 
+: 
+    public fieldsDataBase
+{
+private:
+
+    const   shape&       shape_;
+
+protected:
+    
+    /// check if pointField name exists in Time or time folder 
+    bool pointFieldNameExists(const word& name)const override;
+ 
+
+    /// Loads a pointField with a given name to the Time object. 
+    /// For simulation, it just checks if the name exists 
+    bool loadPointFieldToTime(const word& name) override;
+
+    /// Loads pointStructure to the Time object
+    /// For simulation, it just checks if the name exists
+    bool loadPointStructureToTime() override;
+
+    const shape& getShape() const override;  
+
+    word getPointFieldType(const word& name)const override;
+    
+public:
+
+    TypeInfo("fieldsDataBase<simulation>");
+
+    simulationFieldsDataBase(
+        systemControl& control, 
+        const dictionary& postDict,
+        bool inSimulation,
+        timeValue startTime);
+
+    ~simulationFieldsDataBase() override = default;
+
+    add_vCtor
+    (
+        fieldsDataBase,
+        simulationFieldsDataBase,
+        bool
+    );
+
+    const pointStructure& pStruct()const override;
+
+};
+
+
+}
+
+
+#endif //__simulationFieldsDataBase_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp

@@ -0,0 +1,20 @@
+#include "PostprocessOperationAvMassVelocity.hpp"
+
+pFlow::PostprocessOperationAvMassVelocity::PostprocessOperationAvMassVelocity
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    PostprocessOperationAverage
+    (
+        opDict,
+        opDict.getValOrSet<word>("velocityName", "velocity"),
+        opDict.getValOrSet<word>("massName", "mass"),
+        "all",
+        regPoints,
+        fieldsDB
+    )
+{
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp

@@ -0,0 +1,173 @@
+/*------------------------------- 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 __PostprocessOperationAvMassVelocity_hpp__
+#define __PostprocessOperationAvMassVelocity_hpp__
+
+/*!
+ * @class PostprocessOperationAvMassVelocity
+ * @brief A class for averaging field values within specified regions during post-processing.
+ * 
+ * @details
+ * The PostprocessOperationAvMassVelocity class is a specialized post-processing operation that 
+ * calculates the average of field values within specified regions. It inherits from the 
+ * postprocessOperation base class and implements a weighted averaging operation that 
+ * can be applied to scalar (real), vector (realx3), and tensor (realx4) fields.
+ * 
+ * The average operation follows the mathematical formula:
+ * \f[
+ *    \text{result} = \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * Where:
+ * - \f$ i \f$ represents all particles within the specified processing region
+ * - \f$ j \f$ belongs to a subset of \f$ i \f$ based on an includeMask
+ * - \f$ w_i \f$ is the weight factor for particle \f$ i \f$
+ * - \f$ \phi_i \f$ is the value from the phi field for particle \f$ i \f$
+ * - \f$ \text{field}_j \f$ is the value from the target field for particle \f$ j \f$
+ * 
+ * The calculation can optionally be divided by the region volume (when divideByVolume is set to yes),
+ * which allows calculating normalized averages:
+ * \f[
+ *    \text{result} = \frac{1}{V_{\text{region}}} \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * The averaging can be further filtered using an includeMask to selectively include only
+ * specific particles that satisfy certain criteria.
+ * 
+ * This class supports the following field types:
+ * - real (scalar values)
+ * - realx3 (vector values)
+ * - realx4 (tensor values)
+ * 
+ * @section usage Usage Example
+ * Below is a sample dictionary showing how to configure and use this class:
+ *
+ * ```
+ * processMethod    arithmetic;     // method of performing the sum (arithmetic, uniformDistribution, GaussianDistribution)
+ * processRegion    sphere;         // type of region on which processing is performed 
+ * 
+ * sphereInfo
+ * {
+ *     radius 0.01;
+ *     center (-0.08 -0.08 0.015);
+ * }
+ * 
+ * timeControl default; 
+ * 
+ * /// all the post process operations to be done 
+ * operations
+ * (
+ *     // computes the arithmetic mean of particle velocity 
+ *     averageVel
+ *     {
+ *         function         average;
+ *         field            velocity;
+ *         dividedByVolume  no;     // default is no
+ *         threshold        3;      // default is 1
+ *         includeMask      all;    // include all particles in the calculation
+ *     }
+ *     
+ *     // computes the fraction of par1 in the region 
+ *     par1Fraction
+ *     {
+ *         function        average;
+ *         field           one;     // the "one" field is special - all members have value 1.0
+ *         phi             one;     // default is "one"
+ *         dividedByVolume no;
+ *         includeMask     lessThan;
+ *         
+ *         // diameter of par1 is 0.003, so these settings 
+ *         // will select only particles of type par1
+ *         lessThanInfo
+ *         {
+ *             field      diameter;
+ *             value      0.0031;  
+ *         }
+ *     }
+ * );
+ * ```
+ *
+ * @section defaults Default Behavior
+ * - By default, `phi` is set to the field named "one" which contains value 1.0 for all entries
+ * - `dividedByVolume` is set to "no" by default
+ * - `threshold` is set to 1 by default
+ * - `includeMask` can be set to various filters, with "all" being the default to include all particles
+ *
+ * @section special Special Fields
+ * The field named "one" is a special field where all members have the value 1.0. This makes it
+ * particularly useful for calculating:
+ * 
+ * 1. Volume or number fractions (as shown in the par1Fraction example)
+ * 2. Simple counts when used with an appropriate mask
+ * 3. Normalizing values by particle count
+ *
+ * @see postprocessOperation
+ * @see executeAverageOperation
+ */
+
+#include <variant>
+#include <vector>
+
+#include "PostprocessOperationAverage.hpp"
+#include "regionField.hpp"
+#include "includeMask.hpp"
+
+namespace pFlow
+{
+
+
+class PostprocessOperationAvMassVelocity
+:
+    public PostprocessOperationAverage
+{
+   
+public:
+
+    TypeInfo("PostprocessOperation<avMassVelocity>");
+
+    /// @brief Constructs average operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
+    PostprocessOperationAvMassVelocity(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    /// destructor 
+    ~PostprocessOperationAvMassVelocity() override = default;
+
+    /// add this virtual constructor to the base class
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationAvMassVelocity,
+        dictionary
+    );
+
+}; 
+
+
+}
+
+#endif //__PostprocessOperationAvMassVelocity_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -0,0 +1,140 @@
+#include "PostprocessOperationAverage.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+#include "operationFunctions.hpp"
+
+/// Constructs average processor and initializes result field based on input field type
+pFlow::PostprocessOperationAverage::PostprocessOperationAverage
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(opDict, regPoints, fieldsDB),
+    calculateFluctuation2_(opDict.getValOrSet<Logical>("fluctuation2", Logical(false)))
+{
+    if( fieldType() == getTypeName<real>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, real(0)));
+    }
+    else if( fieldType() == getTypeName<realx3>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx3(0)));
+    }
+    else if( fieldType() == getTypeName<realx4>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx4(0)));
+    }
+    else
+    {
+        fatalErrorInFunction<<" in dictionary "<< opDict.globalName()
+            << " field type is not supported for average operation"
+            << " field type is "<< fieldType()
+            << endl;
+        fatalExit;
+    }
+}
+
+pFlow::PostprocessOperationAverage::PostprocessOperationAverage
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName,
+    const word &includeName, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(opDict, fieldName, phiName, includeName, regPoints, fieldsDB),
+    calculateFluctuation2_(opDict.getValOrSet<Logical>("fluctuation2", Logical(false)))
+{
+    if( fieldType() == getTypeName<real>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, real(0)));
+    }
+    else if( fieldType() == getTypeName<realx3>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx3(0)));
+    }
+    else if( fieldType() == getTypeName<realx4>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx4(0)));
+    }
+    else
+    {
+        fatalErrorInFunction<<" in dictionary "<< opDict.globalName()
+            << " field type is not supported for average operation"
+            << " field type is "<< fieldType()
+            << endl;
+        fatalExit;
+    }
+}
+
+/// Performs weighted average of field values within each region
+bool pFlow::PostprocessOperationAverage::execute
+(
+    const std::vector<span<real>>& weights,
+    const regionField<real>& volFactor
+)
+{
+    auto allField = database().updateFieldAll(fieldName());
+    auto phi = database().updateFieldReal(
+        phiFieldName());
+
+    auto mask = getMask();
+    word procName = processedFieldName();
+    const auto& regP = regPoints();
+    bool dbVol = divideByVolume();
+
+    processedRegFieldPtr_ = makeUnique<processedRegFieldType>
+    (   
+        std::visit([&](auto&& field)->processedRegFieldType
+            { 
+                return executeAverageOperation(
+                procName,
+                field,
+                volFactor,
+                dbVol,
+                weights,
+                phi,
+                mask);
+            },
+            allField)
+    );
+    
+    if(calculateFluctuation2_)
+    {
+        auto& processedRegField = processedRegFieldPtr_();
+        fluctuation2FieldPtr_ = makeUnique<processedRegFieldType>
+        (
+            std::visit([&](auto& field)->processedRegFieldType
+            { 
+                using T = typename std::decay_t<std::remove_reference_t< decltype(field)>>::valueType;
+                if constexpr( std::is_same_v<T,real> ||
+                              std::is_same_v<T,realx3>||
+                              std::is_same_v<T,realx4>)
+                {
+                    return executeFluctuation2Operation(
+                        procName,
+                        field,
+                        std::get<regionField<T>>(processedRegField),
+                        volFactor,
+                        dbVol,
+                        weights,
+                        mask);
+                }
+            },
+            allField)
+        );
+    }
+    
+
+    return true;
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

@@ -0,0 +1,205 @@
+/*------------------------------- 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 __PostprocessOperationAverage_hpp__
+#define __PostprocessOperationAverage_hpp__
+
+/*!
+ * @class PostprocessOperationAverage
+ * @brief A class for averaging field values within specified regions during post-processing.
+ * 
+ * @details
+ * The PostprocessOperationAverage class is a specialized post-processing operation that 
+ * calculates the average of field values within specified regions. It inherits from the 
+ * postprocessOperation base class and implements a weighted averaging operation that 
+ * can be applied to scalar (real), vector (realx3), and tensor (realx4) fields.
+ * 
+ * The average operation follows the mathematical formula:
+ * \f[
+ *    \text{result} = \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * Where:
+ * - \f$ i \f$ represents all particles within the specified processing region
+ * - \f$ j \f$ belongs to a subset of \f$ i \f$ based on an includeMask
+ * - \f$ w_i \f$ is the weight factor for particle \f$ i \f$
+ * - \f$ \phi_i \f$ is the value from the phi field for particle \f$ i \f$
+ * - \f$ \text{field}_j \f$ is the value from the target field for particle \f$ j \f$
+ * 
+ * The calculation can optionally be divided by the region volume (when divideByVolume is set to yes),
+ * which allows calculating normalized averages:
+ * \f[
+ *    \text{result} = \frac{1}{V_{\text{region}}} \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * The averaging can be further filtered using an includeMask to selectively include only
+ * specific particles that satisfy certain criteria.
+ * 
+ * This class supports the following field types:
+ * - real (scalar values)
+ * - realx3 (vector values)
+ * - realx4 (tensor values)
+ * 
+ * @section usage Usage Example
+ * Below is a sample dictionary showing how to configure and use this class:
+ *
+ * ```
+ * processMethod    arithmetic;     // method of performing the sum (arithmetic, uniformDistribution, GaussianDistribution)
+ * processRegion    sphere;         // type of region on which processing is performed 
+ * 
+ * sphereInfo
+ * {
+ *     radius 0.01;
+ *     center (-0.08 -0.08 0.015);
+ * }
+ * 
+ * timeControl default; 
+ * 
+ * /// all the post process operations to be done 
+ * operations
+ * (
+ *     // computes the arithmetic mean of particle velocity 
+ *     averageVel
+ *     {
+ *         function         average;
+ *         field            velocity;
+ *         dividedByVolume  no;     // default is no
+ *         threshold        3;      // default is 1
+ *         includeMask      all;    // include all particles in the calculation
+ *     }
+ *     
+ *     // computes the fraction of par1 in the region 
+ *     par1Fraction
+ *     {
+ *         function        average;
+ *         field           one;     // the "one" field is special - all members have value 1.0
+ *         phi             one;     // default is "one"
+ *         dividedByVolume no;
+ *         includeMask     lessThan;
+ *         
+ *         // diameter of par1 is 0.003, so these settings 
+ *         // will select only particles of type par1
+ *         lessThanInfo
+ *         {
+ *             field      diameter;
+ *             value      0.0031;  
+ *         }
+ *     }
+ * );
+ * ```
+ *
+ * @section defaults Default Behavior
+ * - By default, `phi` is set to the field named "one" which contains value 1.0 for all entries
+ * - `dividedByVolume` is set to "no" by default
+ * - `threshold` is set to 1 by default
+ * - `includeMask` can be set to various filters, with "all" being the default to include all particles
+ *
+ * @section special Special Fields
+ * The field named "one" is a special field where all members have the value 1.0. This makes it
+ * particularly useful for calculating:
+ * 
+ * 1. Volume or number fractions (as shown in the par1Fraction example)
+ * 2. Simple counts when used with an appropriate mask
+ * 3. Normalizing values by particle count
+ *
+ * @see postprocessOperation
+ * @see executeAverageOperation
+ */
+
+#include <variant>
+#include <vector>
+
+#include "postprocessOperation.hpp"
+#include "regionField.hpp"
+#include "includeMask.hpp"
+
+namespace pFlow
+{
+
+
+class PostprocessOperationAverage
+:
+    public postprocessOperation
+{
+private:
+
+    ///< Flag to calculate fluctuation powered by 2     
+    Logical     calculateFluctuation2_;
+
+    /// Result field containing averages for each region (real, realx3, or realx4)
+    uniquePtr<processedRegFieldType>  processedRegFieldPtr_ = nullptr;
+
+    uniquePtr<processedRegFieldType>  fluctuation2FieldPtr_ = nullptr;
+
+public:
+
+    TypeInfo("PostprocessOperation<average>");
+
+    /// @brief Constructs average operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
+    PostprocessOperationAverage(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    PostprocessOperationAverage(
+        const dictionary& opDict,
+        const word& fieldName,
+        const word& phiName,
+        const word& includeName,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    
+    /// destructor 
+    ~PostprocessOperationAverage() override = default;
+
+    /// add this virtual constructor to the base class
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationAverage,
+        dictionary
+    );
+
+    /// @brief Get the processed field containing regional averages
+    /// @return Const reference to average results
+    const processedRegFieldType& processedField()const override
+    {
+        return processedRegFieldPtr_();
+    }
+
+    /// @brief Execute average operation on field values
+    /// @param weights Weight factors for particles
+    /// @return True if successful
+    bool execute(
+        const std::vector<span<real>>& weights,
+        const regionField<real>& volFactor) override;
+
+}; 
+
+
+}
+
+#endif //__PostprocessOperationAverage_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp

@@ -1,8 +1,9 @@
 #include "PostprocessOperationSum.hpp"
 #include "dictionary.hpp"
 #include "fieldsDataBase.hpp"
-#include "fieldFunctions.hpp"
+#include "operationFunctions.hpp"
 
+/// Constructs sum processor and initializes result field based on input field type
 pFlow::PostprocessOperationSum::PostprocessOperationSum
 (
     const dictionary &opDict, 
@@ -37,9 +38,11 @@ pFlow::PostprocessOperationSum::PostprocessOperationSum
     }
 }
 
+/// Performs weighted sum of field values within each region
 bool pFlow::PostprocessOperationSum::execute
 (
-    const std::vector<span<real>>& weights
+    const std::vector<span<real>>& weights,
+    const regionField<real>& volFactor
 )
 {
     auto allField = database().updateFieldAll(fieldName());
@@ -58,7 +61,7 @@ bool pFlow::PostprocessOperationSum::execute
                 return executeSumOperation(
                 procName,
                 field,
-                regP,
+                volFactor,
                 dbVol,
                 weights,
                 phi,

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp

@@ -21,6 +21,107 @@ Licence:
 #ifndef __PostprocessOperationSum_hpp__
 #define __PostprocessOperationSum_hpp__
 
+/*!
+ * @class PostprocessOperationSum
+ * @brief A class for summing field values within specified regions during post-processing.
+ * 
+ * @details
+ * The PostprocessOperationSum class is a specialized post-processing operation that 
+ * calculates the sum of field values within specified regions. It inherits from the 
+ * postprocessOperation base class and implements a weighted summation operation that 
+ * can be applied to scalar (real), vector (realx3), and tensor (realx4) fields.
+ * 
+ * The sum operation follows the mathematical formula:
+ * \f[
+ *    \text{result} = \sum_{i \in \text{processRegion}} w_i \cdot \phi_i \cdot \text{field}_i
+ * \f]
+ * 
+ * Where:
+ * - \f$ i \f$ represents particles within the specified processing region
+ * - \f$ w_i \f$ is the weight factor for particle \f$ i \f$
+ * - \f$ \phi_i \f$ is the value from the phi field for particle \f$ i \f$
+ * - \f$ \text{field}_i \f$ is the value from the      target field for particle \f$ i \f$
+ * 
+ * The calculation can optionally be divided by the region volume (when divideByVolume is set to yes),
+ * which allows calculating density-like quantities:
+ * \f[
+ *    \text{result} = \frac{1}{V_{\text{region}}} \sum_{i \in \text{processRegion}} w_i \cdot \phi_i \cdot \text{field}_i
+ * \f]
+ * 
+ * The summation can be further filtered using an includeMask to selectively include only
+ * specific particles that satisfy certain criteria.
+ * 
+ * This class supports the following field types:
+ * - real (scalar values)
+ * - realx3 (vector values)
+ * - realx4 (tensor values)
+ * 
+ * @section usage Usage
+ * 
+ * To use the PostprocessOperationSum class in a postprocessDataDict file, the following
+ * parameters can be specified:
+ * 
+ * - function: Must be set to "sum" to use this operation
+ * - field: The name of the field to process (e.g., "velocity", "diameter", "one")
+ *   - Special fields like "one" (constant value 1) are also supported
+ *   - Expressions like "cube(diameter)" can be used for mathematical operations
+ * - dividedByVolume: Whether to divide the sum by the region volume (yes/no, default: no)
+ * - includeMask: Optional mask to filter which particles to include in the calculation
+ * 
+ * @section example Example Configuration
+ * 
+ * Here is an example configuration in the postprocessDataDict file:
+ * 
+ * @code
+ * {
+ *     processMethod    arithmetic; 
+ *     processRegion    line;
+ *     
+ *     // the time interval for executing the post-processing 
+ *     // other options: timeStep, default, and settings
+ *     timeControl simulationTime; 
+ *     startTime   1.0;
+ *     endTime     3.0;
+ *     executionInterval 0.1;
+ *      
+ *     // 10 spheres with radius 0.01 along the straight line defined by p1 and p2 
+ *     lineInfo
+ *     {
+ *         p1 (0 0 0);
+ *         p2 (0 0.15 0.15);
+ *         numPoints 10;
+ *         radius    0.01;
+ *     }
+ *     
+ *     operations
+ *     (
+ *         // computes the number density (particles per unit volume)
+ *         numberDensity
+ *         {
+ *             function       sum;
+ *             field          one;        // constant field with value 1.0
+ *             dividedByVolume yes;       // divide by region volume
+ *         }
+ *         
+ *         // computes an approximation of volume fraction
+ *         volumeDensity
+ *         {
+ *             function        sum;
+ *             field           cube(diameter); // d^3, although it differs by pi/6
+ *             dividedByVolume yes;
+ *         }
+ *     );
+ * }
+ * @endcode
+ * 
+ * In this example:
+ * - numberDensity: Calculates the number of particles per unit volume
+ * - volumeDensity: Calculates an approximation of the volume fraction using d³
+ * 
+ * @see postprocessOperation
+ * @see executeSumOperation
+ */
+
 #include <variant>
 #include <vector>
 
@@ -37,21 +138,26 @@ class PostprocessOperationSum
     public postprocessOperation
 {
 private:
-           
+    /// Result field containing sums for each region (real, realx3, or realx4)
-    /// Pointer to the include mask used for masking operations.
     uniquePtr<processedRegFieldType>  processedRegField_ = nullptr;
 
 public:
 
     TypeInfo("PostprocessOperation<sum>");
 
+    /// @brief Constructs sum operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
     PostprocessOperationSum(
         const dictionary& opDict,
         const regionPoints& regPoints,
         fieldsDataBase& fieldsDB);
     
+    /// destructor 
     ~PostprocessOperationSum() override = default;
 
+    /// add this virtual constructor to the base class
     add_vCtor
     (
         postprocessOperation,
@@ -59,16 +165,23 @@ public:
         dictionary
     );
 
+    /// @brief Get the processed field containing regional sums
+    /// @return Const reference to sum results
     const processedRegFieldType& processedField()const override
     {
         return processedRegField_();
     }
 
-    bool execute(const std::vector<span<real>>& weights) override;
+    /// @brief Execute sum operation on field values
+    /// @param weights Weight factors for particles
+    /// @return True if successful
+    bool execute(
+        const std::vector<span<real>>& weights,
+        const regionField<real>& volFactor) override;
 
 }; 
 
 
 }
 
-#endif //__PostprocessOperation_hpp__
+#endif //__PostprocessOperationSum_hpp__

src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp

@@ -18,8 +18,8 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#ifndef __fieldFunctions_hpp__
+#ifndef __operationFunctions_hpp__
-#define __fieldFunctions_hpp__
+#define __operationFunctions_hpp__
 
 #include <vector>
 
@@ -31,27 +31,28 @@ Licence:
 namespace pFlow
 {
 
-
 template<typename T>
 regionField<T> executeSumOperation
 (
     const word&                     regFieldName,
     const span<T>&                  field, 
-    const regionPoints&             regPoints,
+    const regionField<real>&        volFactor,
     const bool                      devideByVol,
     const std::vector<span<real>>&  weights,
     const span<real>&               phi,
     const includeMask::Mask&        mask
 )
 {
+    const auto& regPoints = volFactor.regPoints(); 
     regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)
     {
         auto partIndices = regPoints.indices(reg);
-        auto vols = regPoints.volumes();
+        
         auto w = weights[reg];
-        T sum{};
+        T sum = T{};
         uint n = 0;
         for(auto index:partIndices)
         {
@@ -63,7 +64,7 @@ regionField<T> executeSumOperation
         }
         if(devideByVol)
         {
-            processedField[reg] = sum/vols[reg];
+            processedField[reg] = sum/(volFactor[reg] * vols[reg]);
         }
         else
         {
@@ -80,33 +81,109 @@ regionField<T> executeAverageOperation
 (
     const word&                    regFieldName,
     const span<T>&                 field, 
-    const regionPoints&            regPoints, 
+    const regionField<real>&       volFactor,
+    const bool                     devideByVol,
     const std::vector<span<real>>& weights,
     const span<real>&              phi,
     const includeMask::Mask&       mask
 )
 {
+    
+    const auto& regPoints = volFactor.regPoints(); 
     regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)
     {
         auto partIndices = regPoints.indices(reg);
         auto w = weights[reg];
-        T sumNum{};
+        T sumNum = T{};
-        real sumDen{};
+        real sumDen = 0;
         uint n = 0;
         for(auto index:partIndices)
         {
-            if( index!= -1 && mask( index ))
+            if( index!= -1)
             {
-                sumNum += w[n] * field[index]* phi[index];
+                if( mask(index))
+                {
+                    sumNum += w[n] * field[index]* phi[index];
+                }
+                sumDen += w[n] * phi[index];
             }
-            sumDen += w[n] * phi[index];
+            
             n++;
         }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / (volFactor[reg] * vols[reg]);
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
+        
+    }
 
-        sumDen = max(sumDen, smallValue);
+    return processedField;
-        processedField[reg] = sumNum/sumDen;
+}
+
+
+template<typename T>
+regionField<T> executeFluctuation2Operation
+(
+    const word&                    regFieldName,
+    const span<T>&                 field,
+    const regionField<T>&          fieldAvg,
+    const regionField<real>&       volFactor, 
+    const bool                     devideByVol, 
+    const std::vector<span<real>>& weights,
+    const includeMask::Mask&       mask
+)
+{
+    const auto& regPoints = fieldAvg.regPoints();
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto w = weights[reg];
+        auto vol = volFactor[reg] * vols[reg];
+        T avField{};
+        if(devideByVol)
+        {
+            avField = vol * fieldAvg[reg];
+        }
+        else
+        {
+            avField = fieldAvg[reg];
+        }
+        
+        T sumNum = T{};
+        real sumDen = 0;
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( index!= -1)
+            {
+                if( mask(index))
+                {
+                    sumNum += w[n] * pow( avField- field[index],static_cast<real>(2));
+                }
+                sumDen += w[n];
+            }
+            n++;
+        }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / vol;
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
         
     }
 
@@ -115,4 +192,4 @@ regionField<T> executeAverageOperation
 
 } // namespace pFlow
 
-#endif //__fieldFunctions_hpp__
+#endif //__operationFunctions_hpp__

src/PostprocessData/operation/includeMask/IncludeMask.hpp

@@ -18,6 +18,39 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+/**
+ * @brief A template class implementing includeMask for filtering data based on field values
+ * 
+ * The IncludeMask class creates boolean masks that identify which elements from a field
+ * satisfy a given condition. It applies an operator to each element of a field and
+ * generates a mask (vector of booleans) where true means the element satisfies the
+ * condition and should be included.
+ * 
+ * @tparam T The data type of the field (real, realx3, realx4)
+ * @tparam Operator The operation class that defines the condition (e.g., greaterThanOp)
+ * 
+ * The class has a specialized implementation for allOp operator which includes all elements.
+ * 
+ * Usage example in postprocessDataDict:
+ * ```
+ * // Create a dictionary with the required configuration for filtering
+ * 
+ * {
+ *     includeMask lessThan;
+ *
+ *     // Diameter of par1 is 0.003, so these settings
+ *     // will select only particles of type par1
+ *     lessThanInfo
+ *     {
+ *         field   diameter;
+ *         value   0.0031;
+ *     }
+ *    
+ * }
+ * 
+ * ```
+ */
+
 #ifndef __IncludeMask_hpp__
 #define __IncludeMask_hpp__
 
@@ -31,7 +64,6 @@ Licence:
 namespace pFlow
 {
 
-
 template<typename T, typename Operator>
 class IncludeMask
 :
@@ -39,18 +71,24 @@ class IncludeMask
 {
 public:
 
+    /// Type alias for the mask container returned by getMask()
     using Mask = typename includeMask::Mask;
 
 private:
 	
+    /// Boolean vector the filtering status of each element (true = include)
     std::vector<bool>       mask_;
 
+    /// Comparison operator instance that evaluates the filtering condition
     Operator                operator_;
 
+    /// Name of the field to apply filtering on
     word                    fieldName_;
     
+    /// Timestamp when mask was last updated (-1 indicates never updated)
     timeValue               lastUpdated_ = -1;
 
+    /// Updates the mask based on current field values if needed, returns true if successful
     bool updateMask()
     {
         timeValue t = database().currentTime();
@@ -91,6 +129,7 @@ private:
         return true;
     }
 
+    /// Returns the name of the operator as a string (from operator's TYPENAME)
     static
     word operatorName()
     {
@@ -101,6 +140,7 @@ public:
 
 	TypeInfoTemplate12("IncludeMask", T, Operator);
 
+    /// Constructs an IncludeMask using settings from dictionary and field database
 	IncludeMask(
 		const dictionary& dict, 
         fieldsDataBase& feildsDB)
@@ -113,12 +153,37 @@ public:
                 operatorName()+"Info"
             ).getVal<word>("field"))
 	{}
+    
+    IncludeMask(
+        const word& type,
+		const dictionary& dict, 
+        fieldsDataBase& feildsDB)
+	:
+		includeMask(type, dict, feildsDB),
+		operator_(dict.subDict(operatorName()+"Info")),
+        fieldName_(
+            dict.subDict
+            (
+                operatorName()+"Info"
+            ).getVal<word>("field"))
+	{}
 
-	add_vCtor(
+    /// Add virtual constructor pattern for creating instances
+	add_vCtor
+    (
 		includeMask,
 		IncludeMask,
-		dictionary);
+		dictionary
+    );
 
+    add_vCtor
+    (
+		includeMask,
+		IncludeMask,
+		word
+    );
+
+    /// Returns the mask for filtering elements (updates the mask if necessary)
     Mask getMask() override
     {
         updateMask();
@@ -161,20 +226,40 @@ public:
 
 	TypeInfoTemplate12("IncludeMask", T, allOp<int8>);
 
-	IncludeMask(
+    IncludeMask(
-		const dictionary& opDict, 
+        const dictionary& opDict, 
         fieldsDataBase& feildsDB)
-	:
+    :
-		includeMask(opDict, feildsDB)
+        includeMask(opDict, feildsDB)
-	{
+    {
         span<realx3>  s = database().updatePoints();        
         mask_.resize(s.size(), true);
-	}
+    }
 
-	add_vCtor(
+    IncludeMask(
-		includeMask,
+        const word& type,
-		IncludeMask,
+        const dictionary& opDict, 
-		dictionary);
+        fieldsDataBase& feildsDB)
+    :
+        includeMask(type, opDict, feildsDB)
+    {
+        span<realx3>  s = database().updatePoints();        
+        mask_.resize(s.size(), true);
+    }
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        dictionary
+    );
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        word
+    );
 
     Mask getMask()override
     {

src/PostprocessData/operation/includeMask/includeMask.cpp

@@ -33,11 +33,21 @@ pFlow::includeMask::includeMask
 	database_(fieldDB)
 {}
 
+pFlow::includeMask::includeMask
+(
+	const word &type, 
+	const dictionary &opDict, 
+	fieldsDataBase &fieldsDB
+)
+:
+	database_(fieldsDB)
+{
+}
 
 pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
 (
     const dictionary& opDict, 
-    fieldsDataBase& feildsDB
+    fieldsDataBase& fieldsDB
 )
 {
 	word mask = opDict.getValOrSet<word>("includeMask", "all");
@@ -47,7 +57,7 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
 		auto& maskDict = opDict.subDict(mask+"Info");
 		word maskField = maskDict.getVal<word>("field");
 		
-		if( !feildsDB.getPointFieldType(maskField, fieldType) )
+		if( !fieldsDB.getFieldType(maskField, fieldType) )
 		{
 			fatalErrorInFunction<<"Error in retriving the type of field"
 				<< maskField <<" from dictionary "
@@ -68,7 +78,7 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
 	{
 		auto objPtr = 
 			dictionaryvCtorSelector_[method]
-			(opDict, feildsDB);
+			(opDict, fieldsDB);
 		return objPtr;
 	}
 	else
@@ -87,5 +97,56 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
 	return nullptr;
 }
 
-
+pFlow::uniquePtr<pFlow::includeMask> 
+	pFlow::includeMask::create
+(
+	const word &type, 
+	const dictionary &opDict, 
+	fieldsDataBase &fieldsDB
+)
+{
+	word fieldType;
+	if( type != "all")
+	{
+		auto& maskDict = opDict.subDict(type+"Info");
+		word maskField = maskDict.getVal<word>("field");
+		
+		if( !fieldsDB.getFieldType(maskField, fieldType) )
+		{
+			fatalErrorInFunction<<"Error in retriving the type of field"
+				<< maskField <<" from dictionary "
+				<< maskDict.globalName()
+				<< endl;
+			fatalExit;
+			return nullptr;
+		}
+	}
+	else
+	{
+		fieldType = getTypeName<real>();
+	}
 	
+	word method = angleBracketsNames2("IncludeMask", fieldType, type);
+	
+	if( wordvCtorSelector_.search(method) )
+	{
+		auto objPtr = 
+			wordvCtorSelector_[method]
+			(type, opDict, fieldsDB);
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< 
+			method << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		fatalExit;
+		return nullptr;
+	}
+	return nullptr;
+}

src/PostprocessData/operation/includeMask/includeMask.hpp

@@ -18,6 +18,27 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+/**
+ * @class includeMask
+* @brief Base class for creating inclusion masks for data filtering
+* 
+ * includeMask is an abstract base class for creating masks that filter data
+ * from a fieldsDataBase. Derived classes implement specific filtering criteria
+ * through the getMask() method which returns a Mask object - a wrapper around
+ * a vector of booleans that indicates which elements to include/exclude.
+ * 
+ * This class follows a factory pattern with create() methods that instantiate
+ * the appropriate derived mask type based on dictionary settings.
+ * 
+ * Derived classes can implement various filtering strategies such as:
+ * - Filtering by field values
+ * - Filtering by spatial regions
+ * - Filtering by predefined criteria
+ * 
+ * The Mask objects returned by getMask() can be used in postprocessing operations
+ * to include only the data points that match the specified criteria.
+ */
+
 #ifndef __includeMask_hpp__
 #define __includeMask_hpp__
 
@@ -28,36 +49,49 @@ Licence:
 namespace pFlow
 {
 
-
+// forward declaration
 class fieldsDataBase;
 class dictionary;
 
+
 class includeMask
 {
 public:
 
+	/// @brief Wrapper around a boolean vector that represents elements to include/exclude
+	/// Provides a functional interface to access the underlying mask vector
+	/// where true indicates inclusion and false indicates exclusion.
 	class Mask
 	{
+		/// @brief Reference to the underlying boolean vector
 		const std::vector<bool>& mask_;
 
 		public:
 
+		/// @brief Constructor from a boolean vector
+		/// @param msk Boolean vector where true means include, false means exclude
 		Mask(const std::vector<bool>& msk)
 		:
 			mask_(msk)
 		{}
 
+		/// @brief Copy constructor
 		Mask(const Mask&) = default;
 
+		/// @brief Move constructor
 		Mask(Mask&&) = default;
 				
+		/// @brief Destructor
 		~Mask()=default;
 
+		/// @brief Get the size of the mask
 		auto size()const
 		{
 			return mask_.size();
 		}
 
+		/// @brief Check if element at index i should be included
+		/// @param i Index to check
 		bool operator()(uint32 i)const
 		{
 			return mask_[i];
@@ -66,16 +100,28 @@ public:
 
 private:
 	
+	/// @brief Reference to the database containing fields to be filtered
     fieldsDataBase&     database_;
 	
 public:
 
 	TypeInfo("includeMask");
 
+	/// @brief Constructor
+	/// @param opDict Dictionary containing operation settings
+	/// @param feildsDB Database of fields 
 	includeMask(const dictionary& opDict, fieldsDataBase& feildsDB);
 
+	/// @brief Constructor with explicit type
+	/// @param type Type of mask to create
+	/// @param opDict Dictionary containing operation settings
+	/// @param feildsDB Database of fields 
+    includeMask(const word& type, const dictionary& opDict, fieldsDataBase& feildsDB);
+
+	/// @brief Virtual destructor
 	virtual ~includeMask() = default;
 
+	/// @brief Virtual constructor pattern implementation using dictionary
     create_vCtor
     (
         includeMask,
@@ -85,24 +131,58 @@ public:
         ),
         (opDict, feildsDB)		
     );
+
+	/// @brief Virtual constructor pattern implementation using type and dictionary
+    create_vCtor
+    (
+        includeMask, 
+        word,
+        (
+            const word& type, 
+            const dictionary& opDict, 
+            fieldsDataBase& feildsDB
+        ),
+        (type, opDict, feildsDB)
+    );
 	 
+	/// @brief Get const access to the database
+	/// @return Const reference to the fields database
 	const fieldsDataBase& database()const
 	{
 		return database_;
 	}
 
+	/// @brief Get non-const access to the database
+	/// @return Reference to the fields database
     fieldsDataBase& database()
     {
         return database_;
     }
 
+	/// @brief Get the mask for filtering elements
+	/// @return Mask object indicating which elements to include
     virtual 
     Mask getMask()= 0;
 
+	/// @brief Factory method to create appropriate mask type from dictionary
+	/// @param opDict Dictionary with mask settings including type
+	/// @param feildsDB Database of fields to filter
+	/// @return Unique pointer to created mask instance
 	static 
 	uniquePtr<includeMask> create(
         const dictionary& opDict, 
         fieldsDataBase& feildsDB);
+    
+	/// @brief Factory method to create mask with explicit type
+	/// @param type Type of mask to create
+	/// @param opDict Dictionary with mask settings
+	/// @param feildsDB Database of fields to filter
+	/// @return Unique pointer to created mask instance
+    static
+    uniquePtr<includeMask> create(
+        const word& type,
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB);
 	
 };
 

src/PostprocessData/operation/includeMask/maskOperations.hpp

@@ -1,3 +1,22 @@
+/*------------------------------- 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 __maskOperation_hpp__
 #define __maskOperation_hpp__

src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp

@@ -82,14 +82,33 @@ bool writeField
 
 }
 
-
-
 pFlow::postprocessOperation::postprocessOperation
 (
     const dictionary &opDict, 
     const regionPoints& regPoints,
     fieldsDataBase &fieldsDB
 )
+:
+    postprocessOperation
+    (
+        opDict,
+        opDict.getVal<word>("field"),
+        opDict.getValOrSet<word>("phi", "one"),
+        opDict.getValOrSet<word>("includeMask", "all"),
+        regPoints,
+        fieldsDB
+    )
+{}
+
+pFlow::postprocessOperation::postprocessOperation
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName,
+    const word& includeName,  
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
 :
     operationDict_(opDict),
     threshold_
@@ -98,7 +117,7 @@ pFlow::postprocessOperation::postprocessOperation
     ),
     divideByVolume_
     (
-        opDict.getValOrSet<Logical>("dividedByVolume", Logical(false))
+        opDict.getValOrSet<Logical>("divideByVolume", Logical(false))
     ),
     regionPoints_
     (
@@ -110,25 +129,24 @@ pFlow::postprocessOperation::postprocessOperation
     ),
     fieldName_
     (
-        opDict.getValOrSet<word>("field", "one")
+        fieldName
     ),
     phiFieldName_
     (
-        opDict.getValOrSet<word>("phi", "one")
+        phiName
     ),
     includeMask_
     (
-        includeMask::create(opDict, fieldsDB)
+        includeMask::create(includeName, opDict, fieldsDB)
     )
 {
 
-    if(!fieldsDB.getPointFieldType(fieldName_, fieldType_))
+    if(!fieldsDB.getFieldType(fieldName_, fieldType_))
     {
         fatalErrorInFunction;
         fatalExit;
     }
 }
-
 const pFlow::Time& pFlow::postprocessOperation::time() const
 {
     return database_.time();
@@ -141,7 +159,7 @@ bool pFlow::postprocessOperation::write(const fileSystem &parDir) const
     if(!osPtr_)
     {
         fileSystem path = parDir+(
-            processedFieldName() + ".Start_" + ti.prevTimeName());
+            processedFieldName() + ".Start_" + ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
         
         regPoints().write(osPtr_());

src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp

@@ -20,6 +20,53 @@ Licence:
 #ifndef __postprocessOperation_hpp__
 #define __postprocessOperation_hpp__
 
+/*! 
+ *  @class postprocessOperation
+ *  @file postprocessOperation.hpp
+ *  @brief Base class for post-processing operations on particle data.
+ *         This class provides the foundational structure and functionality
+ *         for performing various post-processing operations on simulation data.
+ *
+ *  @details
+ *  The postprocessOperation class operates on field data (specified in the input dictionary)
+ *  and performs specific operations on that field within defined regions. It serves as
+ *  part of the post-processing framework in phasicFlow to analyze particle simulation results.
+ *  
+ *  Operations are performed on specific subsets of particles defined by region points and
+ *  can be filtered using include masks. The class supports different field types (real,
+ *  realx3, realx4) through the processedRegFieldType variant.
+ *  
+ *  The main operations supported include:
+ *
+ *  1. Sum operation:
+ *     - Calculates: 
+ *       \f[
+ *          \text{result} = \sum_{i \in \text{processRegion}} w_i \cdot \phi_i \cdot \text{field}_i
+ *       \f]
+ *     - Where \f$ i \f$ belongs to the particles in the specified processRegion
+ *     - \f$ w_i \f$ is the weight factor for particle \f$ i \f$
+ *     - \f$ \phi_i \f$ is the value from the phi field for particle \f$ i \f$
+ *     - \f$ \text{field}_i \f$ is the value from the target field for particle \f$ i \f$
+ *     - Implemented in the derived class PostprocessOperationSum
+ *
+ *  2. Average operation:
+ *     - Calculates: 
+ *       \f[
+ *          \text{result} = \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                              {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ *       \f]
+ *     - Where \f$ i \f$ belongs to all particles in the specified processRegion
+ *     - \f$ j \f$ belongs to a subset of \f$ i \f$ based on an includeMask defined in the input dictionary
+ *     - This allows calculating regional averages on specific subsets of particles
+ *
+ *  The class uses threshold values to exclude regions with insufficient particles
+ *  and supports optional division by volume for density-like calculations. Results are written
+ *  to files for later analysis or visualization.
+ *
+ *  @note The actual processing is performed by derived classes that implement
+ *  the execute() method for specific operation types.
+ */
+
 #include <variant>
 
 #include "virtualConstructor.hpp"
@@ -33,6 +80,9 @@ Licence:
 namespace pFlow
 {
 
+/// Type alias for processed region field types.
+/// Only regionField<real>, regionField<realx3>, and regionField<realx4> are supported
+/// in the postprocessOperation class.
 using processedRegFieldType = std::variant
     <
         regionField<real>,
@@ -40,14 +90,10 @@ using processedRegFieldType = std::variant
         regionField<realx4>
     >;
 
+/// - forward declaration
 class fieldsDataBase;
 class Time;
 
-/*!
- *  @brief Base class for post-processing operations.
- *         This class provides the basic structure and functionality
- *         for performing post-processing operations on simulation data.
- */
 class postprocessOperation
 {
 public:
@@ -88,16 +134,31 @@ private:
     
 public:
 
-
+    /// Type info
     TypeInfo("postprocessOperation");
 
+    /// Constructor
+    /// @param opDict Dictionary containing operation-specific parameters.
+    /// @param regPoints Reference to the region points used in the operation.
+    /// @param fieldsDB Reference to the fields database containing field data.
     postprocessOperation(
         const dictionary& opDict,
         const regionPoints& regPoints,
         fieldsDataBase& fieldsDB );
-
+    
+    postprocessOperation(
+        const dictionary& opDict,
+        const word& fieldName,
+        const word& phiName,
+        const word& includeName,  
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB
+    );
+    
+    /// destructor
     virtual ~postprocessOperation()=default;
-
+    
+    /// Active the virtual constructor for creating derived classes.
     create_vCtor(
         postprocessOperation,
         dictionary,
@@ -108,74 +169,102 @@ public:
         ),
         (opDict, regPoints, fieldsDB));
     
+    /// Access to regionPoints instance 
     const regionPoints& regPoints()const
     {
         return regionPoints_;
     }
 
+    /// Access to fields database instance
     const fieldsDataBase& database()const
     {
         return database_;
     }
 
+    /// Access to fields database instance
     fieldsDataBase& database()
     {
         return database_;
     }
 
+    /// Access to the time instance
     const Time& time()const;
 
+    /// Return the name of the processed field.
     word processedFieldName()const
     {
         return operationDict_.name();
     }
     
+    /// return the name of the field to be processed.
     const word& fieldName()const
     {
         return fieldName_;
     }
 
+    /// return the type name of the field to be processed.
     const word& fieldType()const
     {
         return fieldType_;
     }
 
+    /// return the name of the phi field to be processed.
     const word& phiFieldName()const
     {
         return phiFieldName_;
     }
 
+    /// Access to the operation dictionary 
     const dictionary& operationDict()const
     {
         return operationDict_;
     }
+
+    /// return threshold value 
+    /// which is used to exclude the regions which contain
+    /// particles fewer than this value.
     const uint32 threshold()const
     {
         return threshold_;
     }
 
+    /// whether the result is divided by volume of the region 
     bool divideByVolume()const
     {
         return divideByVolume_();
     }
 
+    /// return the include mask
     Mask getMask()
     {
         return includeMask_().getMask();
     }
 
+    /// return the processed field
     virtual
     const processedRegFieldType& processedField()const=0;
 
-    virtual 
+    /// execute the operation
-    bool execute(const std::vector<span<real>>& weights) = 0;
+    /// @param weights Vector of weights for the operation.
+    /// @param volFactor a factor to be multiplied by the volume of the region
+    virtual bool execute(
+        const std::vector<span<real>>& weights, 
+        const regionField<real>& volFactor) = 0;
 
+    /// write the result to a file 
+    /// @param parDir Parent directory for the output file.
     virtual 
     bool write(const fileSystem &parDir)const;
 
+    /// 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;}
 
+    /// Create the polymorphic object using the virtual constructor.
+    /// @param opDict Dictionary containing operation-specific parameters.
+    /// @param regPoints Reference to the region points used in the operation.
+    /// @param fieldsDB Reference to the fields database containing field data.
     static
     uniquePtr<postprocessOperation> create(
         const dictionary& opDict,

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -36,6 +36,12 @@ pFlow::PostprocessComponent<RegionType,ProcessMethodType>::PostprocessComponent
     (
         regionPointsPtr_().size()
     ),
+    volumeFactor_
+    (
+        "volumeFactor",
+        regionPointsPtr_(),
+        1.0
+    ),
     operationDicts_(readDictList("operations", dict))
 {
 
@@ -67,6 +73,10 @@ bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::execute
         return true;
     }
 
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
     // update processing methods
     auto& regPoints = this->regPoints();
     
@@ -102,7 +112,7 @@ bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::execute
     
     for(auto& op:operatios_)
     {
-        if( !op->execute(weights) )
+        if( !op->execute(weights, volumeFactor_) )
         {
             fatalErrorInFunction
                 <<"error occured in executing operatoin defined in dict "

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp

@@ -51,6 +51,8 @@ private:
     /// Method for processing the selected particles data 
     std::vector<ProcessMethodType>  regionsProcessMethod_;
 
+    regionField<real>               volumeFactor_;
+
     bool    executed_{false};
 
     dictionaryList                  operationDicts_;
@@ -62,6 +64,16 @@ protected:
         return regionsProcessMethod_;
     }
 
+    regionField<real>& volumeFactor()
+    {
+        return volumeFactor_;
+    }
+    
+    const regionField<real>& volumeFactor()const
+    {
+        return volumeFactor_;
+    }
+
 public:
 
     // type info 

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp

@@ -51,10 +51,12 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
             regs[i] = arithmetic(); // Changed from uniformDistribution() to arithmetic()
+            volFactor[i] = 1.0;
         }
     }
 

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp

@@ -23,6 +23,7 @@ Licence:
 
 #include "PostprocessComponent.hpp"
 #include "GaussianDistribution.hpp"
+#include "numericConstants.hpp"
 
 namespace pFlow
 {
@@ -51,15 +52,20 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
+
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
-            regs[i] = GaussianDistribution(c[i], pow(d[i],2));
+            auto r = d[i]/2;
+            regs[i] = GaussianDistribution(c[i], pow(r/3.0,2));
+            volFactor[i] = 0.677683 / (4.0/3.0*Pi*r);
         }
     }
 
     // add the virtual constructor
-    add_vCtor(
+    add_vCtor
+    (
         postprocessComponent,
         PostprocessComponentGaussian,
         dictionary

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp

@@ -51,10 +51,12 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
             regs[i] = uniformDistribution();
+            volFactor[i] = 1.0;
         }
     }
 

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp

@@ -24,10 +24,19 @@ Licence:
 
 // region types 
 #include "sphereRegionPoints.hpp"
-
+#include "lineRegionPoints.hpp"
+#include "multipleSpheresRegionPoints.hpp"
 
 
 template class pFlow::PostprocessComponentGaussian<pFlow::sphereRegionPoints>;
 template class pFlow::PostprocessComponentUniform<pFlow::sphereRegionPoints>;
 template class pFlow::PostprocessComponentArithmetic<pFlow::sphereRegionPoints>;
 
+template class pFlow::PostprocessComponentGaussian<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::multipleSpheresRegionPoints>;
+
+template class pFlow::PostprocessComponentGaussian<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::lineRegionPoints>;
+

src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp

@@ -92,7 +92,12 @@ bool pFlow::particleProbePostprocessComponent::execute
         executed_ = false;
         return true;
     }
-    
+
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
+
     if(!regionPointsPtr_().update())
     {
         fatalErrorInFunction
@@ -130,7 +135,7 @@ bool pFlow::particleProbePostprocessComponent::write(const fileSystem& parDir)co
     if( !osPtr_)
     {
         // file is not open yet
-        fileSystem path = parDir + (name_+".Start_"+ti.prevTimeName());
+        fileSystem path = parDir + (name_+".Start_"+ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
         regionPointsPtr_().write(osPtr_());
     }

src/PostprocessData/postprocessData/postprocessData.cpp

@@ -22,15 +22,18 @@ Licence:
 #include "List.hpp"
 #include "systemControl.hpp"
 #include "postprocessData.hpp"
-#include "fileDictionary.hpp"
 #include "postprocessGlobals.hpp"
 #include "postprocessComponent.hpp"
 
-pFlow::postprocessData::postprocessData(const systemControl &control)
+pFlow::postprocessData::postprocessData
+(
+    const systemControl &control,
+    timeValue startTime
+)
 :
     auxFunctions(control),
+    inSimulation_(startTime<0.0?true:false),
     time_(control.time()),
-    fieldsDataBase_(control.time(), true),
     dict_
     (
         objectFile
@@ -40,16 +43,27 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
             objectFile::READ_IF_PRESENT,
             objectFile::WRITE_NEVER
         )
-    ),
+    )
-    componentsDicts_(readDictList("components", dict_))
 {
     postProcessGlobals::defaultDir__ = CWD()/pFlow::postProcessGlobals::defaultRelDir__;
     
     // if dictionary is not provided, no extra action is required. 
-    if( !dict_.fileExist() )
+    if( !dict_.fileExist() || !dict_.headerOk() )
     {
+        WARNING<<"You requested postprocessData function while,"
+            <<" the dictionary system/postprocessDataDict does not exist."
+            <<" This feature is disabled in the current run."<<END_WARNING;
         return;
     }
+    
+    fieldsDataBasePtr_= fieldsDataBase::create
+        (
+            const_cast<systemControl&>(control),
+            dict_,
+            inSimulation_,
+            startTime
+        );
+    
 
     activeInSimulation_ = dict_.getValOrSet<Logical>(
         "activeInSimulation", 
@@ -72,11 +86,13 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
             "execution");
     }
     
-    for(auto& compDict:componentsDicts_)
+    componentsDictsPtr_ = makeUnique<dictionaryList>(readDictList("components", dict_));
+
+    for(auto& compDict:*componentsDictsPtr_)
     {
         postprocesses_.push_back( postprocessComponent::create(
             compDict, 
-            fieldsDataBase_, 
+            fieldsDataBasePtr_(), 
             defaultTimeControlPtr_() ));
     }
 
@@ -84,14 +100,16 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
 
 bool pFlow::postprocessData::execute() 
 {
+    if( inSimulation_ && !activeInSimulation_() ) return true; 
+
     const auto& ti = time_.TimeInfo();
 
     for(auto& component:postprocesses_)
     {
-        if(!component->execute(ti))
+        if(!component->execute(ti, !inSimulation_) )
         {
             fatalErrorInFunction
-                <<"Error occured in executing postprocess component: "
+                <<"Error occurred in executing postprocess component: "
                 <<component->name()<<endl;
             return false;
         }
@@ -102,12 +120,15 @@ bool pFlow::postprocessData::execute()
 
 bool pFlow::postprocessData::write() const
 {
+    if( inSimulation_ && !activeInSimulation_() ) return true; 
+    
     for(auto& component:postprocesses_)
     {   
         if(!component->executed())
         {
             continue;
         }
+
         if(!component->write(postProcessGlobals::defaultDir__/component->name()))
         {
             fatalErrorInFunction
@@ -118,3 +139,8 @@ bool pFlow::postprocessData::write() const
     }
     return true;
 }
+
+void pFlow::postprocessData::setOutputDirectory(const fileSystem &path) const
+{
+    postProcessGlobals::defaultDir__ = path;
+}

src/PostprocessData/postprocessData/postprocessData.hpp

@@ -20,14 +20,14 @@ Licence:
 #ifndef __postprocessData_hpp__
 #define __postprocessData_hpp__
 
-#include "auxFunctions.hpp"
 #include "Logical.hpp"
 #include "ListPtr.hpp"
 #include "fileDictionary.hpp"
-#include "baseTimeControl.hpp"
+#include "dictionaryList.hpp"
+#include "auxFunctions.hpp"
 #include "fieldsDataBase.hpp"
 #include "postprocessComponent.hpp"
-#include "dictionaryList.hpp"
+
 
 namespace pFlow 
 {
@@ -36,6 +36,7 @@ class systemControl;
 class Time;
 class timeInfo;
 
+
 /**
  * @class postprocessData
  * @brief An interface class for handling post-processing of simulation data.
@@ -47,7 +48,11 @@ class postprocessData
 :
   public auxFunctions
 {
-    /// Indicates if a post-processing is active during simulatoin 
+    /// Indicates if this is post-processing during simulation or 
+    /// post-simulation 
+    bool                    inSimulation_;
+
+    /// Indicates if a post-processing is active during simulation 
     Logical                 activeInSimulation_{false};
 
     /// a list of active post-process components 
@@ -57,13 +62,13 @@ class postprocessData
     const Time&           time_;
 
     /// Database for all the points fields on the host 
-    fieldsDataBase        fieldsDataBase_;
+    uniquePtr<fieldsDataBase>        fieldsDataBasePtr_;
 
     /// file dictionary that is constructed from the file (postProcessDataDict)
     fileDictionary        dict_;
 
     /// list of dictionaries for postprocess components 
-    dictionaryList        componentsDicts_;
+    uniquePtr<dictionaryList>        componentsDictsPtr_ = nullptr;
 
     /// @brief  default time control that can be used for all post-process components
     uniquePtr<baseTimeControl> defaultTimeControlPtr_= nullptr;
@@ -76,7 +81,7 @@ public:
     ///  this constructor is used when postprocesing is active
     ///  during simulation.
     /// @param control const reference to systemControl 
-    postprocessData(const systemControl& control);
+    postprocessData(const systemControl& control, timeValue startTime = -1.0);
 
     ~postprocessData()override = default;
 
@@ -89,11 +94,19 @@ public:
 
     bool execute() override;
     
-    
-     
     bool write()const override;
     
+    fieldsDataBase& database()
+    {
+        return fieldsDataBasePtr_();
+    }
 
+    const fieldsDataBase& database()const
+    {
+        return fieldsDataBasePtr_();
+    }
+
+    void setOutputDirectory(const fileSystem& path)const;
 };
 
 } // namespace pFlow

src/PostprocessData/postprocessData/postprocessGlobals.hpp

@@ -27,7 +27,7 @@ namespace pFlow::postProcessGlobals
 {
 
 static fileSystem defaultDir__;
-inline const word defaultRelDir__ = "VTK/postprocessData";
+inline const word defaultRelDir__ = "postprocessData";
 
 }