Update readme.md

corrections for readme.md file postprocessing

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;
-objectType     dicrionary;
+objectName     shapes;
+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;    

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;
@@ -69,13 +85,15 @@ public:
 
 	realVector boundingDiameter()const override;
 	
-		real coarseGrainFactor(uint32 index)const ;
+	realVector volume()const override;
+
+	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

@@ -64,7 +64,8 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		"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;
@@ -75,3 +106,36 @@ 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

@@ -61,8 +61,27 @@ public:
 		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
+(
+    systemControl& control,
+    const dictionary& postDict,
+    bool inSimulation,
+    timeValue startTime
+)
 :
-    time_(time), 
+    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,
-    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);
-        word space;
-        if(!pointFieldGetType(fieldTypeName, originalType, space))
+        // the name is in the Time object 
+        if(pointFieldNameExists(pointFieldName))
         {
-            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, outputType))
+    if(!inputOutputType(func, originalType, typeAfterFunction))
     {
-        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, 
-    word &originalType, 
-    word &typeAfterFunction
-)
+    const word& compoundName,
+    word& originalType, 
+    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)
 {
-    
-    bool shouldUpdate = checkForUpdate("position", forceUpdate);
+    const word fName = "position";
+    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>(
-            "position", 
-            pstruct.activePointsHost());
+        allFields_.emplaceBackOrReplace<PointsTypeHost>
+        (
+            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 
-            static_cast<const pointStructure&>(
-                time_.lookupObject<dynamicPointStructure>(pointStructureFile__)
-            );
+        dbType = "fieldsDataBase<simulation>";
     }
     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 <concepts>
-
+#include "fieldsDataBaseDclr.hpp"
 
 #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:
 
@@ -100,123 +72,111 @@ private:
 
     // - Member variables
 
-        /// const reference to the Time object
-        const Time&         time_;
-
         /// List to store all the point fields 
         anyList             allFields_;
 
         /// Map to store the last capture time of each field
         wordMap<timeValue>  captureTime_;
 
-        /// 
+        /// Reference to the Time object
+        Time&               time_; 
+
+        /// Flag indicating if we're in simulation mode
         bool                inSimulation_ = false;
 
-        static const inline std::map<word, word> reservedFieldNames_
+        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
         {
-            {"position", "realx3"},
-            {"one", "real"}
-        };
+            return shapeType_;
+        }
 
-    // - Private member functions
+        /// 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);
 
-        /// @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 
+        template<ValidFieldType T>
+        span<T> updateReservedField(const word& name, bool forceUpdate = false);
+        
         span<real> createOrGetRealField(const word& name);
 
-        /**
-         * @brief Parses a compound field name to extract the base field name and function.
-         * 
-         * 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
-         * name (e.g., "velocity") and the function to be applied (e.g., `Functions::Magnitude`).
-         * 
-         * The function supports the following syntax for compound field names:
-         * - `fieldName` (no function applied)
-         * - `functionName(fieldName)`
-         * 
-         * Supported function names are defined in the `Functions` enum.
-         * 
-         * @param compoundFieldName The compound field name to parse.
-         * @param fieldName A reference to a `word` where the extracted base field name
-         *                  will be stored.
-         * @param func A reference to a `Functions` enum where the identified function
-         *             will be stored. If no function is applied, this will be set to
-         *             `Functions::None`.
-         * 
-         * @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.
-         */
+        span<real> createOrGetVolume(bool forceUpdate=false);
+
+        span<real> createOrGetDensity(bool forceUpdate=false);
+
+        span<real> createOrGetOne(bool forceUpdate=false);
+
+        span<real> createOrGetMass(bool forceUpdate=false);
+
+        span<real> createOrGetI(bool forceUpdate=false);
+
+        /// Map of reserved field names with their corresponding types
+        static const inline std::map<word, word> reservedFieldNames_
+        {
+            {"position", "realx3"},
+            {"one", "real"},
+            {"volume", "real"},
+            {"density", "real"},
+            {"mass", "real"},
+            {"I", "real"}
+        };
+        
         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,7 +191,21 @@ 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
@@ -243,128 +217,49 @@ public:
             return time_;
         }
         
+        Time& time()
+        {
+            return time_;
+        }
+
     // - Public Member Functions
 
-        /**
-         * @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(
+        bool getFieldTypeNameFunction
+        (
             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(
+        
+        bool getFieldType
+        (
             const word& compoundName,
             word& originalType, 
-            word& typeAfterFunction);
+            word& typeAfterFunction
+        ) const;
 
-        /// overload for function getPointFieldType without `originalType` argument
-        bool getPointFieldType(
+        bool getFieldType
+        (
             const word& compoundName,
-            word& typeAfterFunction);
+            word& typeAfterFunction
+        ) const;
 
-        /**
-         * @brief Updates the points data and returns a span to the updated points.
-         * 
-         * This function ensures that the points data is up-to-date by checking if an update
-         * is necessary. If the data is outdated or if a forced update is requested, it retrieves
-         * the latest points data and stores it in the internal fields database. The function
-         * then returns a span to the updated points data for further use.
-         * 
-         * @param forceUpdate A boolean flag indicating whether to force an update of the points
-         *                    data regardless of its current state. Defaults to `false`.
-         * 
-         * @return A span of `realx3` representing the updated points data.
-         */
+        
+        /// update pointStructure if necessary 
         span<realx3> updatePoints(bool forceUpdate = false);
         
-        /**
-         * @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
-         */
+        /// update a field with realx3 type 
         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);
@@ -379,10 +274,44 @@ public:
             const word& compoundName, 
             bool forceUpdate = false);
         
-        const pointStructure& pStruct()const;
+        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;
+        }
+
+        /// 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 updatePoints(forceUpdate);            
-            }
-            else
-            {
-                fatalErrorInFunction<< "Error in getting the type name of field: "<<
-                    name<<", with type name: "<< getTypeName<T>() <<endl;
-                fatalExit;
-                return span<T>(nullptr, 0);
-            }
+            return updateReservedField<T>(name, true);
         }   
         else
         {
-            const auto& pField = time_.lookupObject<pointField_D<T>>(name);
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
-            name, 
-            pField.activeValuesHost());
+            if( loadPointFieldToTime(name) )
+            {
+                const auto& pField = time_.lookupObject<pointField_D<T>>(name);
+                allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
+                name, 
+                pField.activeValuesHost());
+            }
+            else
+            {
+                fatalErrorInFunction<<"Error in loading the pointField "<<name<<endl;
+                fatalExit;
+            }
+            
         }
     }
 
@@ -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:
-           
-    /// Pointer to the include mask used for masking operations.
+    /// Result field containing sums for each region (real, realx3, or realx4)
     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__
-#define __fieldFunctions_hpp__
+#ifndef __operationFunctions_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{};
-        real sumDen{};
+        T sumNum = T{};
+        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++;
         }
         
-        sumDen = max(sumDen, smallValue);
-        processedField[reg] = sumNum/sumDen;
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / (volFactor[reg] * vols[reg]);
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
+        
+    }
+
+    return processedField;
+}
+
+
+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)
@@ -114,11 +154,36 @@ public:
             ).getVal<word>("field"))
 	{}
     
-	add_vCtor(
+    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 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(
-		const dictionary& opDict, 
+    IncludeMask(
+        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,
-		dictionary);
+    IncludeMask(
+        const word& type,
+        const dictionary& opDict, 
+        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,
@@ -86,24 +132,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 
-    bool execute(const std::vector<span<real>>& weights) = 0;
+    /// execute the operation
+    /// @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

@@ -93,6 +93,11 @@ bool pFlow::particleProbePostprocessComponent::execute
         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,17 +43,28 @@ 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", 
         Logical{true});
@@ -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";
 
 }
 

src/PostprocessData/processMethod/GaussianDistribution.hpp

@@ -26,7 +26,7 @@ Licence:
 #include "typeInfo.hpp"
 #include "types.hpp"
 #include "span.hpp"
-
+#include "numericConstants.hpp"
 
 namespace pFlow
 {
@@ -74,15 +74,15 @@ public:
         for(uint32 i=0; i<indices.size(); i++)
         {
             auto x = points[indices[i]]-center;
-            auto f = exp(- dot(x,x)/(2*variance_));
+            auto f = exp(- dot(x,x)/(2*variance_))/sqrt(2.0*Pi*variance_);
             weight_[i] = f;
             sum += f;
         }
 
-        for(auto& w: weight_)
+        /*for(auto& w: weight_)
         {
             w /= sum;
-        }   
+        } */  
         return true;
     }
 

src/PostprocessData/readme.md

@@ -0,0 +1,363 @@
+# PostprocessData Module in phasicFlow
+
+The `PostprocessData` module in phasicFlow provides powerful tools for analyzing particle-based simulations both during runtime (in-simulation) and after simulation completion (post-simulation). This document explains how to configure and use the postprocessing features through the dictionary-based input system.
+
+## Overview
+
+Postprocessing in phasicFlow allows you to:
+
+- Extract information about particles in specific regions of the domain
+- Calculate statistical properties such as averages and sums of particle attributes
+- Track specific particles throughout the simulation
+- Apply different weighing methods when calculating statistics
+- Perform postprocessing at specific time intervals
+
+## Setting Up Postprocessing
+
+Postprocessing is configured through a dictionary file named `postprocessDataDict` which should be placed in the `settings` directory. Below is a detailed explanation of the configuration options.
+
+### Basic Configuration
+
+The input dictionary, **settings/postprocessDataDict**, may look like this:
+
+```cpp
+// PostprocessData dictionary
+
+
+// Enable/disable postprocessing during simulation
+runTimeActive yes;  // Options: yes, no
+
+// Shape type - only needed when doing post-simulation processing
+shapeType sphere;   // Options depend on the simulation type: sphere, grain, etc.
+
+// Default time control for postprocessing components
+defaultTimeControl
+{
+    timeControl         timeStep;  // Options: timeStep, simulationTime, settings
+    startTime           0;         // Start time for postprocessing
+    endTime             1000;      // End time for postprocessing
+    executionInterval   150;       // How frequently to run postprocessing
+}
+
+// List of postprocessing components
+components
+(
+    // Component definitions here...
+);
+```
+
+If you want to activate in-simulaiton postprocessing, you need to add these lines to the `settings/settingsDict` file:
+
+```cpp
+libs   ("libPostprocessData.so");
+
+auxFunctions    postprocessData;
+```
+
+This will link the postprocessing library to your simulation, allowing you to use its features. Note that, anytime you want to deactivate the in-simulation postprocessing, you can simply change the `runTimeActive` option to `no` in `postprocessDataDict` file.
+
+## Time Control Options
+
+Each postprocessing component can either use the default time control settings or define its own. There are three main options for time control:
+
+| Option | Description | Required Parameters |
+|--------|-------------|---------------------|
+| `timeStep` | Controls execution based on simulation time steps | `startTime`, `endTime`, `executionInterval` |
+| `simulationTime` | Controls execution based on simulation time | `startTime`, `endTime`, `executionInterval` |
+| `settings` | Uses parameters from settingsDict file | None (defined elsewhere) |
+| `default` | Uses the default time control settings (uses `defaultTimeControl` settings)| None (uses default) |
+
+If no time control is specified, the `default` option is used automatically.
+
+## Processing Methods
+
+The postprocessing module provides several methods for processing particle data. They are categorized into two main groups: bulk and individual methods.
+
+- **Bulk Methods**: Operate on all particles that are located in a specified locations/regions (cells, spheres, etc.).
+- **Individual Methods**: Operate on specific particles, allowing for targeted particle property extraction.
+
+| Method | Property type | Description | Formula |
+|--------|------------------|-------------|---------|
+| `arithmetic` | bulk | Simple arithmetic mean/sum with equal weights | Each particle contributes equally |
+| `uniformDistribution` | bulk | Each particle contributes inversely proportional to the total number of particles | $w_i = 1/n$ where $n$ is the number of particles |
+| `GaussianDistribution` | bulk | Weight contribution based on distance from center with Gaussian falloff | $w_i = \exp(-\|x_i - c\|^2/(2\sigma^2))/\sqrt{2\pi\sigma^2}$ |
+| `particleProbe` | individual | Extracts values from specific particles | Direct access to particle properties |
+
+## Region Types
+
+Regions define where in the domain the postprocessing operations are applied:
+
+| Region Type | Description | Required Parameters | Compatible with |
+|-------------|-------------|---------------------|-----------------|
+| `sphere` | A spherical region | `radius`, `center` | bulk |
+| `multipleSpheres` | Multiple spherical regions | `centers`, `radii` | bulk |
+| `line` | Spheres along a line with specified radius | `p1`, `p2`, `nSpheres`, `radius` | bulk |
+| `centerPoints` | Specific particles selected by ID | `ids` | individual |
+
+## Processing Operations
+
+Within each processing region of type `bulk`, you can define multiple operations to be performed:
+
+### Available Functions
+
+| Function | Property type | Description | Formula | Required Parameters |
+|----------|---------------|-------------|---------|---------------------|
+| `average` | bulk | Weighted average of particle field values | see Equation 1 | `field`, `phi` (optional), `threshold` (optional), `includeMask` (optional), `divideByVolume` (optional) |
+| `sum` | bulk | Weighted sum of particle field values | see Equation 2 | `field`, `phi` (optional),`threshold` (optional), `includeMask` (optional), `divideByVolume` (optional) |
+
+Equation 1:
+
+$$\text{result} = \frac{\sum_j w_j \cdot \phi_j \cdot \text{field}_j}{\sum_i w_i \cdot \phi_i}$$
+
+Equation 2:
+
+$$\text{result} = \sum_j w_j \cdot \phi_j \cdot \text{field}_j$$
+
+where:
+
+- $j$ is the index of the particle in the region that also satisfies the `includeMask`
+- $i$ is the index of the particle in the region
+- $w_j$ is the weight of particle $j$ based on the selected processing method
+- $\phi_j$ is the value of the `phi` field for particle $j$ (default is 1)
+- $field_j$ is the value of the specified field for particle $j$
+
+
+### Derived Functions
+
+In addition to the above basic functions, some derived functions are available for specific calculations:
+
+| Function | Property type | Description | Formula | Required Parameters |
+|----------|---------------|-------------|---------|---------------------|
+|`avMassVelocity` | bulk | Average velocity weighted by mass | $\frac{\sum_{i \in \text{region}} m_i \cdot v_i}{\sum_{i \in \text{region}} m_i}$ | - |
+
+### Available Fields
+
+All the pointFields in the simulation database (for in-simulation processing), or the ones stored in the time folders (for post-simulation processing) can be referenced in the operations. In addition to them, some extra fields are available for use in the operations. The following fields are available for use in the operations:
+
+1. Extra fileds to be used in post-processing operations:
+
+| Field | Field Type | Description | Default Value |
+|-------|------------|-------------|---------------|
+| `position` | `realx3` | Particle positions | - |
+| `one` | `real` | Value 1 for each particle | 1 |
+| `mass` | `real` | Particle mass | - |
+| `density` | `real` | Particle density | - |
+| `volume` | `real` | Particle volume | - |
+| `diameter` | `real` | Particle diameter | - |
+| `I` | `real` | Moment of inertia | - |
+
+2. Common fields which are available in the simulation database/time folders:
+
+| Field | Field Type | Description |
+|-------|------------|-------------|
+| `velocity` | `realx3` | Particle velocity |
+| `rVelocity` | `realx3` | Particle rotational velocity |
+| `acceleration` | `realx3` | Particle acceleration |
+| `rAcceleration` | `realx3` | Particle rotational acceleration |
+| `contactForce` | `realx3` | Particle contact force |
+| `contactTorque` | `realx3` | Particle contact torque |
+| `id` | `integer` | Particle ID |
+| `shapeIndex` | `integer` | Particle shape index |
+
+The above fields may vary from one type of simulation to other. Pleas note that this is only a tentative list.
+
+### Optional Parameters
+
+| Parameter | Description | Default | Options |
+|-----------|-------------|---------|---------|
+| `divideByVolume` | Divide result by region volume | `no` | `yes`, `no` |
+| `threshold` | Exclude regions with fewer particles | 1 | Integer value |
+| `includeMask` | Filter particles based on a field value | `all` | `all`, `lessThan`, `greaterThan`, `between`, `lessThanOrEq`, `greaterThanEq`, `betweenEq` |
+
+## Examples
+
+### Example 1: Probing Individual Particles
+
+```cpp
+velocityProb
+{
+    processMethod    particleProbe;
+    processRegion    centerPoints;
+    selector         id;
+    field            component(position,y);
+    ids              (0 10 100);
+    timeControl      default;
+}
+```
+
+This example extracts the y-component of the position for particles with IDs 0, 10, and 100.
+
+### Example 2: Processing in a Spherical Region
+
+```cpp
+on_single_sphere
+{
+    processMethod    arithmetic;
+    processRegion    sphere;
+    
+    sphereInfo
+    {
+        radius  0.01;
+        center  (-0.08 -0.08 0.015);
+    }
+    
+    timeControl default;
+    
+    operations
+    (
+        averageVel
+        {
+            function        average;
+            field           mag(velocity);
+            divideByVolume  no;
+            threshold       3;
+            includeMask     all;
+        }
+        
+        par1Fraction
+        {
+            function        average;
+            field           one;
+            phi             one;
+            divideByVolume  no;
+            includeMask     lessThan;
+            
+            lessThanInfo
+            {
+                field       diameter;
+                value       0.0031;
+            }
+        }
+        
+        numberDensity
+        {
+            function        sum;
+            field           one;
+            phi             one;
+            divideByVolume  yes;
+        }
+    );
+}
+```
+
+This example defines a sphere region and performs three operations:
+
+1. Calculate the average of velocity magnitude of particles
+2. Calculate the fraction of particles with diameter less than 0.0031
+3. Calculate the number density by summing and dividing by volume
+
+### Example 3: Processing Along a Line
+
+In this example, a line region is defined. The `lineInfo` section specifies the start and end points of the line, the number of spheres to create along the line, and the radius of each point. Bulk properties are calculated in each sphere, based on the properties of particles contained in each sphere. 
+
+```cpp
+along_a_line
+{
+      processMethod    arithmetic;
+      processRegion    line;
+      
+      timeControl simulationTime;
+      startTime   1.0;
+      endTime     3.0;
+      executionInterval 0.1;
+      
+      lineInfo
+      {
+        p1   (0 0 0);
+        p2   (0 0.15 0.15);
+        nSpheres    10;
+        radius    0.01;
+      }
+      
+      operations
+      (
+          bulkDensity
+          {
+              function       sum;
+              field          mass;
+              phi            one;
+              divideByVolume yes;
+          }
+          
+          volumeDensity
+          {
+            function         sum;
+            field            volume;
+            divideByVolume   yes;
+          }
+     );
+}
+```
+
+This example creates 10 spherical regions along a line from (0,0,0) to (0,0.15,0.15) and calculates the bulk density and volume density in each region.
+
+## Advanced Features
+
+### Special functions applied on fields
+
+You can access specific components of vector fields (`realx3`) using the `component` function:
+
+```cpp
+field  component(position,y);  // Access the y-component of position
+```
+
+Here is a complete list of these special functions: 
+
+| Function Name | Valid field Type | Example |
+|-----------------|------------|---------|
+| `component` | `realx3` | `component(velocity,x)` |
+| `abs` | `real` | `abs(s)` |
+| `square` | `real` | `square(I)` |
+| `cube` | `real` | `cube(diameter)` |
+| `squre root` | `real` | `sqrt(mass)` |
+| `magnitude` | `realx3` | `mag(contactForce)` |
+| `magnitude square` | `realx3` | `magSquare(velocity)` |
+| `magnitude cube` | `realx3` | `magCube(velocity)` |
+| `magnitude square root` | `realx3` | `magSqrt(acceleration)` |
+
+
+### Particle Filtering with includeMask
+
+The `includeMask` parameter allows you to filter particles based on field values:
+
+```cpp
+includeMask  lessThan;
+lessThanInfo
+{
+    field   diameter;
+    value   0.0031;
+}
+```
+
+Supported masks:
+
+- `all`: Include all particles (default)
+- `lessThan`: Include particles where field < value
+- `greaterThan`: Include particles where field > value
+- `between`: Include particles where value1 < field < value2
+- `lessThanOrEq`: Include particles where field ≤ value
+- `greaterThanOrEq`: Include particles where field ≥ value
+- `betweenEq`: Include particles where value1 ≤ field ≤ value2
+
+## Implementation Notes
+
+- The postprocessing system can work both during simulation (`runTimeActive yes`) or after simulation completion.
+- When using post-simulation mode, you must specify the correct `shapeType` to properly initialize the shape objects.
+- Results are written to output files in the case directory with timestamps.
+- The `threshold` parameter helps eliminate statistical noise in regions with few particles.
+- Setting `divideByVolume` to `yes` normalizes results by the volume of the region, useful for calculating densities.
+
+## Mathematical Formulations
+
+For weighted `bulk` properties calculation:
+
+$$ \text{average} = \frac{\sum_{i \in \text{region and includeMask}} w_i \cdot \phi_i \cdot \text{field}_i}{\sum_{i \in \text{region}} w_i \cdot \phi_i} $$
+
+For weighted summing:
+
+$$ \text{sum} = \sum_{i \in \text{region and includeMask}} w_i \cdot \phi_i \cdot \text{field}_i $$
+
+If `divideByVolume` is set to `yes`, the result is divided by the volume of the region:
+
+$$ \text{volumetric result} = \frac{\text{result}}{V_{\text{region}}} $$
+