Merge pull request #222 from PhasicFlow/main

from main

.github/scripts/sync-wiki.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python3
+
+import os
+import re
+import yaml
+import sys
+
+# Constants
+REPO_URL = "https://github.com/PhasicFlow/phasicFlow"
+REPO_PATH = os.path.join(os.environ.get("GITHUB_WORKSPACE", ""), "repo")
+WIKI_PATH = os.path.join(os.environ.get("GITHUB_WORKSPACE", ""), "wiki")
+MAPPING_FILE = os.path.join(REPO_PATH, ".github/workflows/markdownList.yml")
+
+def load_mapping():
+    """Load the markdown to wiki page mapping file."""
+    try:
+        with open(MAPPING_FILE, 'r') as f:
+            data = yaml.safe_load(f)
+            return data.get('mappings', [])
+    except Exception as e:
+        print(f"Error loading mapping file: {e}")
+        return []
+
+def convert_relative_links(content, source_path):
+    """Convert relative links in markdown content to absolute URLs."""
+    # Find markdown links with regex pattern [text](url)
+    md_pattern = r'\[([^\]]+)\]\(([^)]+)\)'
+    
+    # Find HTML img tags
+    img_pattern = r'<img\s+src=[\'"]([^\'"]+)[\'"]'
+    
+    def replace_link(match):
+        link_text = match.group(1)
+        link_url = match.group(2)
+        
+        # Skip if already absolute URL or anchor
+        if link_url.startswith(('http://', 'https://', '#', 'mailto:')):
+            return match.group(0)
+        
+        # Get the directory of the source file
+        source_dir = os.path.dirname(source_path)
+        
+        # Create absolute path from repository root
+        if link_url.startswith('/'):
+            # If link starts with /, it's already relative to repo root
+            abs_path = link_url
+        else:
+            # Otherwise, it's relative to the file location
+            abs_path = os.path.normpath(os.path.join(source_dir, link_url))
+            if not abs_path.startswith('/'):
+                abs_path = '/' + abs_path
+        
+        # Convert to GitHub URL
+        github_url = f"{REPO_URL}/blob/main{abs_path}"
+        return f"[{link_text}]({github_url})"
+    
+    def replace_img_src(match):
+        img_src = match.group(1)
+        
+        # Skip if already absolute URL
+        if img_src.startswith(('http://', 'https://')):
+            return match.group(0)
+        
+        # Get the directory of the source file
+        source_dir = os.path.dirname(source_path)
+        
+        # Create absolute path from repository root
+        if img_src.startswith('/'):
+            # If link starts with /, it's already relative to repo root
+            abs_path = img_src
+        else:
+            # Otherwise, it's relative to the file location
+            abs_path = os.path.normpath(os.path.join(source_dir, img_src))
+            if not abs_path.startswith('/'):
+                abs_path = '/' + abs_path
+        
+        # Convert to GitHub URL (use raw URL for images)
+        github_url = f"{REPO_URL}/raw/main{abs_path}"
+        return f'<img src="{github_url}"'
+    
+    # Replace all markdown links
+    content = re.sub(md_pattern, replace_link, content)
+    
+    # Replace all img src tags
+    content = re.sub(img_pattern, replace_img_src, content)
+    
+    return content
+
+def process_file(source_file, target_wiki_page):
+    """Process a markdown file and copy its contents to a wiki page."""
+    source_path = os.path.join(REPO_PATH, source_file)
+    target_path = os.path.join(WIKI_PATH, f"{target_wiki_page}.md")
+    
+    print(f"Processing {source_path} -> {target_path}")
+    
+    try:
+        # Check if source exists
+        if not os.path.exists(source_path):
+            print(f"Source file not found: {source_path}")
+            return False
+        
+        # Read source content
+        with open(source_path, 'r') as f:
+            content = f.read()
+        
+        # Convert relative links
+        content = convert_relative_links(content, source_file)
+        
+        # Write to wiki page
+        with open(target_path, 'w') as f:
+            f.write(content)
+        
+        return True
+    
+    except Exception as e:
+        print(f"Error processing {source_file}: {e}")
+        return False
+
+def main():
+    # Check if wiki directory exists
+    if not os.path.exists(WIKI_PATH):
+        print(f"Wiki path not found: {WIKI_PATH}")
+        sys.exit(1)
+    
+    # Load mapping
+    mappings = load_mapping()
+    if not mappings:
+        print("No mappings found in the mapping file")
+        sys.exit(1)
+    
+    print(f"Found {len(mappings)} mappings to process")
+    
+    # Process each mapping
+    success_count = 0
+    for mapping in mappings:
+        source = mapping.get('source')
+        target = mapping.get('target')
+        
+        if not source or not target:
+            print(f"Invalid mapping: {mapping}")
+            continue
+        
+        if process_file(source, target):
+            success_count += 1
+    
+    print(f"Successfully processed {success_count} of {len(mappings)} files")
+    
+    # Exit with error if any file failed
+    if success_count < len(mappings):
+        sys.exit(1)
+
+if __name__ == "__main__":
+    main()

.github/workflows/markdownList.yml

@@ -0,0 +1,18 @@
+# This file maps source markdown files to their target wiki pages
+# format:
+# - source: path/to/markdown/file.md
+#   target: Wiki-Page-Name
+mappings:
+  - source: benchmarks/readme.md
+    target: Performance-of-phasicFlow
+  - source: benchmarks/helicalMixer/readme.md
+    target: Helical-Mixer-Benchmark
+  - source: benchmarks/rotatingDrum/readme.md
+    target: Rotating-Drum-Benchmark
+  - source: doc/mdDocs/howToBuild-V1.0.md
+    target: How-to-build-PhasicFlow‐v‐1.0
+  - source: tutorials/README.md
+    target: Tutorials
+  - source: doc/mdDocs/phasicFlowFeatures.md
+    target: Features-of-PhasicFlow
+  # Add more mappings as needed

.github/workflows/sync-wiki.yml

@@ -0,0 +1,60 @@
+name: Sync-Wiki
+
+on:
+  push:
+    branches:
+      - main
+    paths:
+      - "**/*.md"
+      - ".github/workflows/sync-wiki.yml"
+      - ".github/workflows/markdownList.yml"
+      - ".github/scripts/sync-wiki.py"
+  workflow_dispatch:
+
+jobs:
+  sync-wiki:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Repository
+        uses: actions/checkout@v3
+        with:
+          path: repo
+      - name: Checkout Wiki
+        uses: actions/checkout@v3
+        with:
+          repository: ${{ github.repository }}.wiki
+          path: wiki
+        continue-on-error: true
+      - name: Create Wiki Directory if Not Exists
+        run: |
+          if [ ! -d "wiki" ]; then
+            mkdir -p wiki
+            cd wiki
+            git init
+            git config user.name "${{ github.actor }}"
+            git config user.email "${{ github.actor }}@users.noreply.github.com"
+            git remote add origin "https://github.com/${{ github.repository }}.wiki.git"
+          fi
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: '3.10'
+      - name: Install dependencies
+        run: pip install pyyaml
+      - name: Sync markdown files to Wiki
+        run: |
+          python $GITHUB_WORKSPACE/repo/.github/scripts/sync-wiki.py
+        env:
+          GITHUB_REPOSITORY: ${{ github.repository }}
+      - name: Push changes to wiki
+        run: |
+          cd wiki
+          git config user.name "${{ github.actor }}"
+          git config user.email "${{ github.actor }}@users.noreply.github.com"
+          git add .
+          if git status --porcelain | grep .; then
+            git commit -m "Auto sync wiki from main repository"
+            git push --set-upstream https://${{ github.actor }}:${{ github.token }}@github.com/${{ github.repository }}.wiki.git master -f
+          else
+            echo "No changes to commit"
+          fi

benchmarks/helicalMixer/readme.md

@@ -0,0 +1 @@
+# Helical Mixer Benchmark (phasicFlow v-1.0)

benchmarks/readme.md

@@ -0,0 +1,7 @@
+
+# Benchmarks
+
+Benchmakrs has been done on two different simulations: a simulation with simple geometry (rotating drum) and a simulation with complex geometry (helical mixer).
+
+- [rotating drum](./rotatingDrum/readme.md)
+- [helical mixer](./helicalMixer/readme.md)

benchmarks/rotatingDrum/readme.md

@@ -0,0 +1,96 @@
+# Rotating Drum Benchmark (phasicFlow v-1.0)
+
+## Overview
+
+This benchmark compares the performance of phasicFlow with a well-stablished commercial DEM software for simulating a rotating drum with varying particle counts (250k to 8M particles). The benchmark measures both computational efficiency and memory usage across different hardware configurations.
+
+## Simulation Setup
+
+<div align="center">
+    <img src="./images/commericalDEMsnapshot.png"/>
+    <div align="center">  
+        <p>Figure 1. Commercial DEM simulation snapshot</p>
+    </div>
+</div>
+
+<div align="center">
+    <img src="./images/phasicFlow_snapshot.png"/>
+    <div align="center">  
+        <p>Figure 2. phasicFlow simulation snapshot and visualized using Paraview</p>
+    </div>
+</div>
+
+### Hardware Specifications
+
+<div align="center">
+    Table 1. Hardware specifications used for benchmarking.
+</div>
+
+|  System     |           CPU            |             GPU              | Operating System |
+| :---------: | :----------------------: | :--------------------------: | :--------------: |
+|   Laptop    | Intel i9-13900HX 2.2 GHz | NVIDIA GeForce RTX 4050Ti 6G | Windows 11 24H2  |
+| Workstation | Intel Xeon 4210 2.2 GHz  |     NVIDIA RTX A4000 16G     |   Ubuntu 22.04   |
+
+### Simulation Parameters
+
+<div align="center">
+    Table 2. Parameters for rotating drum simulations.
+</div>
+
+| Case     | Particle Diameter | Particle Count | Drum Length | Drum Radius |
+| :-------: | :---------------: | :--------------: | :------------------: | :------------------: |
+| 250k     | 6 mm              | 250,000        | 0.8 m       | 0.2 m       |
+| 500k     | 5 mm              | 500,000        | 0.8 m       | 0.2 m       |
+| 1M       | 4 mm              | 1,000,000      | 0.8 m       | 0.2 m       |
+| 2M       | 3 mm              | 2,000,000      | 1.2 m       | 0.2 m       |
+| 4M       | 3 mm              | 4,000,000      | 1.6 m       | 0.2 m       |
+| 8M       | 2 mm              | 8,000,000      | 1.6 m       | 0.2 m       |
+
+The time step for all simulations was set to 1.0e-5 seconds and the simulation ran for 4 seconds.
+
+## Performance Comparison
+
+### Execution Time
+
+<div align="center">
+    Table 3. Total calculation time (minutes) for different configurations.
+</div>
+
+|     Software      | 250k   | 500k   | 1M     | 2M     | 4M     | 8M      |
+| :---------------: | :----: | :-----: | :-----: | :-----: | :-----: | :------: |
+| phasicFlow-4050Ti | 54 min | 111 min | 216 min | 432 min | -       | -        |
+| Commercial DEM-4050Ti | 68 min | 136 min | 275 min | 570 min | -       | -        |
+| phasicFlow-A4000  | 38 min | 73 min  | 146 min | 293 min | 589 min | 1188 min |
+
+The execution time scales linearly with particle count. phasicFlow demonstrates approximately:
+
+- 20% faster calculation than the well-established commercial DEM software on the same hardware
+- 30% performance improvement when using the NVIDIA RTX A4000 compared to the RTX 4050Ti
+
+<div align="center">
+    <img src="./images/performance1.png"/>
+    <p>Figure 3. Calculation time comparison between phasicFlow and the well-established commercial DEM software.</p>
+</div>
+
+### Memory Usage
+
+<div align="center">
+    Table 4. Memory consumption for different configurations.
+</div>
+
+|     Software      | 250k   | 500k   | 1M      | 2M      | 4M      | 8M      |
+| :---------------: | :-----: | :-----: | :-----: | :-----: | :-----: | :-----: |
+| phasicFlow-4050Ti | 252 MB  | 412 MB  | 710 MB  | 1292 MB | -       | -       |
+| Commercial DEM-4050Ti | 485 MB | 897 MB | 1525 MB | 2724 MB | -       | -       |
+| phasicFlow-A4000  | 344 MB  | 480 MB  | 802 MB  | 1386 MB | 2590 MB | 4966 MB |
+
+Memory efficiency comparison:
+
+- phasicFlow uses approximately 0.7 GB of memory per million particles
+- Commercial DEM software uses approximately 1.2 GB of memory per million particles
+- phasicFlow shows ~42% lower memory consumption compared to the commercial alternative
+- The memory usage scales linearly with particle count in both software packages. But due to memory limitations on GPUs, it is possible to run larger simulation on GPUs with phasicFlow.
+
+## Run Your Own Benchmarks
+
+The simulation case setup files are available in this folder for users interested in performing similar benchmarks on their own hardware. These files can be used to reproduce the tests and compare performance across different systems.

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

@@ -2,12 +2,14 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName particleInsertion;
-objectType dicrionary;
-fileFormat  ASCII;
+
+objectName     shapes;
+objectType     dictionary;
+fileFormat     ASCII;
 /*---------------------------------------------------------------------------*/
-active         No;                               // is insertion active -> Yes or No
 
-collisionCheck No;                               // is checked          -> Yes or No
+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  60;                    // 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
+    {
+        distance   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

@@ -2,9 +2,14 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */ 
-objectName particleInsertion;
-objectType dicrionary;
-fileFormat  ASCII;
-/*---------------------------------------------------------------------------*/
-active         No;                                      // is checked          -> Yes or No
 
+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  60;                    // 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
+    {
+        distance   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  60;                    // 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
+    {
+        distance   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

@@ -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.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/geometryDict

@@ -35,7 +35,7 @@ surfaces
         
         radius2     0.2;                   // radius at p2
         
-        resolution  24;                    // number of divisions
+        resolution  60;                    // number of divisions
         
         material    wallMat;               // material name of this wall
         

benchmarks/rotatingDrum/rotatingDrum_4mParticles/settings/particlesDict

@@ -27,9 +27,9 @@ positionParticles
 
     orderedInfo
     {
-        diameter   0.003;           // minimum space between centers of particles
+        distance   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  60;                    // 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
+    {
+        distance   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  60;                    // 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
+    {
+        distance   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;    

doc/mdDocs/howToBuild-V1.0.md

@@ -0,0 +1,136 @@
+# How to build PhasicFlow-v-1.0
+
+You can build PhasicFlow for CPU or GPU. You can have a single build or oven multiple builds on a machine. Here you learn how to have a single build of PhasicFlow, in various modes of execution. You can install PhasicFlow-v-1.0 on **Ubuntu-22.04 LTS** and **Ubuntu-24.04 LTS**. Installing it on older versions of Ubuntu needs some additional steps to meet the requirements which are not covered here. 
+
+If you want to install PhasicFlow on **Windows OS**, just see [this page](https://www.cemf.ir/installing-phasicflow-v-1-0-on-ubuntu/) for more information. 
+
+# Required packages
+
+You need a list of packages installed on your computer before building PhasicFlow:
+
+* git, for cloning the code and package management
+* g++, for compiling the code
+* cmake, for generating build system
+* Cuda-12.x or above (if GPU is targeted), for compiling the code for CUDA execution.
+ 
+### Installing packages 
+
+Execute the following commands to install the required packages (Except Cuda). tbb is installed automatically.
+
+```bash
+sudo apt update
+sudo apt install -y git g++ cmake cmake-qt-gui
+```
+
+### Installing Cuda for GPU execution 
+
+If you want to build PhasicFlow to be executed on an nvidia-GPU, you need to install the latest version of Cuda compiler (Version 12.x or above), which is compatible with your hardware and OS, on your computer.
+
+# How to build? 
+
+Here you will learn how to build PhasicFlow for single execution mode. Follow the steps below to install it on your computer. 
+Tested operating systems are:
+* Ubuntu-22.04 LTS
+* Ubuntu-24.04 LTS
+
+### Step 1: Package check
+Make sure that you have installed all the required packages on your computer. See above for more information.
+
+
+### Step 2: Cloning PhasicFlow
+Create the PhasicFlow folder in your home folder and then clone the source code into that folder:
+```bash
+cd ~
+mkdir PhasicFlow
+cd PhasicFlow
+git clone https://github.com/PhasicFlow/phasicFlow.git
+mv phasicFlow phasicFlow-v-1.0
+```
+### Step 3: Environmental variables
+Opne the bashrc file using the following command:
+
+```bash
+$ gedit ~/.bashrc
+```
+
+and add the following line to the end of the file, **save** and **close** it.
+
+```bash
+source $HOME/PhasicFlow/phasicFlow-v-1.0/cmake/bashrc
+```
+
+this will introduce a new source file for setting the environmental variables of PhasicFlow. If you want to load these variables in the current open terminal, you need to source it. Or, simply **close the terminal** and **open a new terminal**.
+
+### Step 4: Building PhasicFlow
+Follow one of the followings to build PhasicFlow for one mode of execution.
+
+#### Serial build for CPU
+In a **new terminal** enter the following commands:
+```bash
+cd ~/PhasicFlow/phasicFlow-v-1.0
+mkdir build
+cd build
+cmake ../ -DpFlow_Build_Serial=On -DCMAKE_BUILD_TYPE=Release
+make install -j4
+```
+For faster builds, use `make install -j`. This will use all the CPU cores on your computer for building. 
+
+#### OpenMP build for CPU
+
+```bash
+cd ~/PhasicFlow/phasicFlow-v-1.0
+mkdir build
+cd build
+cmake ../ -DpFlow_Build_OpenMP=On -DCMAKE_BUILD_TYPE=Release
+make install -j4
+```
+
+#### GPU build for parallel execution on CUDA-enabled GPUs
+
+```bash
+cd ~/PhasicFlow/phasicFlow-v-1.0
+mkdir build
+cd build
+cmake ../ -DpFlow_Build_Cuda=On -DCMAKE_BUILD_TYPE=Release
+cmake ../ -DpFlow_Build_Cuda=On -DCMAKE_BUILD_TYPE=Release
+make install -j4
+```
+
+After building, `bin`, `include`, and `lib` folders will be created in `~/PhasicFlow/phasicFlow-v-1.0/` folder. Now you are ready to use PhasicFlow.
+
+**note 1**: When compiling the code in parallel, you need to have enough RAM on your computer. As a rule, you need 1 GB free RAM per each processor on your computer for compiling in parallel.
+You may want to use fewer number of cores on your computer by using the following command:
+
+```bash
+make install -j3
+```
+
+the above command only uses 3 cores for compiling. 
+
+**note 2**: By default PhasicFlow is compiled with **double** as floating point variable. You can compile it with **float**. Just in the command line of camke added `-DpFlow_Build_Double=Off` flag to compile it with float. For example if you are building for cuda, you can enter the following command:
+
+```bash
+cmake ../ -DpFlow_Build_Cuda=On -DpFlow_Build_Double=Off
+```
+
+### Step 5: Testing
+
+In the current terminal or a new terminal enter the following command:
+
+```bash
+checkPhasicFlow
+```
+
+This command shows the host and device environments and software version. If PhasicFlow was build correctly, you would get the following output:
+
+```
+Initializing host/device execution spaces . . . 
+  Host execution space is Serial
+  Device execution space is Serial
+
+  You are using phasicFlow v-1.0 (copyright(C): www.cemf.ir)
+  In this build, double is used for floating point operations and uint32for indexing.
+  This is not a build for MPI execution
+
+Finalizing host/device execution space ....
+```

doc/mdDocs/howToBuild.md

@@ -1,151 +0,0 @@
-# How to build PhasicFlow {#howToBuildPhasicFlow}
-
-You can build PhasicFlow for CPU or GPU. You can have a single build or oven multiple builds on a machine. Here you learn how to have a single build of PhasicFlow, in various modes of execution. 
-
-# Required packages
-You need a list of packaged installed on your computer before building PhasicFlow:
-* git, for cloning the code and package management
-* g++, for compiling the code
-* cmake, for generating build system
-* tbb, a parallel library for STL algorithms
-* Cuda (if GPU is targeted), for compiling the code for CUDA execution.
-* Kokkos, the parallelization backend of PhasicFlow 
-
-### git 
-if git is not installed on your computer, enter the following commands 
-```
-$ sudo apt update
-$ sudo apt install git
-```
- 
-
-### g++ (C++ compiler)
-The code is tested with g++ (gnu C++ compiler). The default version of g++ on Ubuntu 18.04 LTS or upper is sufficient for compiling. If it is not installed on your operating system, enter the following command:
-
-``` 
-$ sudo apt update
-$ sudo apt install g++
-```
-
-
-### CMake
-You also need to have CMake-3.22 or higher installed on your computer.
-``` 
-$ sudo apt update
-$ sudo apt install cmake
-```
-
-### tbb (2020.1-2 or higher)
-For **Ubuntu 20.04 LTS or higher versions**, you can install tbb using apt. For now, some parallel algorithms on host side rely on tbb parallel library (C++ parallel backend). Use e following commands to install it:
-```
-$ sudo apt update
-$ sudo apt install libtbb-dev
-```
-If you are compiling on **Ubuntu-18.04 LTS**, you need to enter the following commands to get the right version (2020.1-2 or higher) of tbb:
-```
-$ wget "http://archive.ubuntu.com/ubuntu/pool/universe/t/tbb/libtbb2_2020.1-2_amd64.deb"
-$ sudo dpkg --install libtbb2_2020.1-2_amd64.deb
-$ wget "http://archive.ubuntu.com/ubuntu/pool/universe/t/tbb/libtbb-dev_2020.1-2_amd64.deb"
-$ sudo dpkg --install libtbb-dev_2020.1-2_amd64.deb
-```
-### Cuda
-If you want to build PhasicFlow to be executed on an nvidia-GPU, you need to install the latest version of Cuda compiler, which is compatible with your hardware and OS, on your computer. 
-
-# How to build? 
-Here you will learn how to build PhasicFlow for single execution mode. Follow the steps below to install it on your computer. 
-Tested operating systems are:
-* Ubuntu 18.04 LTS
-* Ubuntu 20.04 LTS
-* Ubuntu 22.04 LTS
-
-### Step 1: Package check
-Make sure that you have installed all the required packages on your computer. See above for more information.
-
-### Step 2: Cloning Kokkos
-It is assumed that Kokkos source is located in the home folder of your computer. Clone the latest version of Kokkos into your home folder:
-```
-$ cd ~
-$ mkdir Kokkos
-$ cd Kokkos
-$ git clone https://github.com/kokkos/kokkos.git
-```
-or simply download and extract the source code of Kokkos in `~/Kokkos` folder. In the end, the top level CMakeLists.txt file should be located in `~/Kokkos/kokkos` folder. 
-
-### Step 3: Cloning PhasicFlow
-Create the PhasicFlow folder in your home folder and then clone the source code into that folder:
-```
-$ cd ~
-$ mkdir PhasicFlow
-$ cd PhasicFlow
-$ git clone https://github.com/PhasicFlow/phasicFlow.git
-```
-### Step 4: Environmental variables
-Opne the bashrc file using the following command:
-
-`$ gedit ~/.bashrc`
-
-and add the following line to the end of the file, **save** and **close** it.
-
-`source $HOME/PhasicFlow/phasicFlow/cmake/bashrc`
-
-this will introduce a new source file for setting the environmental variables of PhasicFlow. If you want to load these variables in the current open terminal, you need to source it. Or, simply **close the terminal** and **open a new terminal**.
-
-### Step 5: Building PhasicFlow
-Follow one of the followings to build PhasicFlow for one mode of execution.
-#### Serial build for CPU
-In a **new terminal** enter the following commands:
-```
-$ cd ~/PhasicFlow/phasicFlow
-$ mkdir build
-$ cd build
-$ cmake ../ -DpFlow_Build_Serial=On
-$ make install
-```
-For faster builds, use `make install -j`. This will use all the CPU cores on your computer for building. 
-#### OpenMP build for CPU
-```
-$ cd ~/PhasicFlow/phasicFlow
-$ mkdir build
-$ cd build
-$ cmake ../ -DpFlow_Build_OpenMP=On
-$ make install
-```
-#### GPU build for parallel execution on CUDA-enabled GPUs
-```
-$ cd ~/PhasicFlow/phasicFlow
-$ mkdir build
-$ cd build
-$ cmake ../ -DpFlow_Build_Cuda=On
-$ make install
-```
-
-After building, `bin`, `include`, and `lib` folders will be created in `~/PhasicFlow/phasicFlow/` folder. Now you are ready to use PhasicFlow.
-
-**note 1**: When compiling the code in parallel, you need to have enough RAM on your computer. As a rule, you need 1 GB free RAM per each processor in your computer for compiling in parallel.
-You may want to use fewer number of cores on your computer by using the following command:
-
-`$ make install -j 3`
-
-the above command only uses 3 cores for compiling. 
-
-**note 2**: By default PhasicFlow is compiled with **double** as floating point variable. You can compile it with **float**. Just in the command line of camke added `-DpFlow_Build_Double=Off` flag to compile it with float. For example if you are building for cuda, you can enter the following command:
-
-`$ cmake ../ -DpFlow_Build_Cuda=On -DpFlow_Build_Double=Off`
-
-### Step 6: Testing
-In the current terminal or a new terminal enter the following command:
-
-`$ checkPhasicFlow`
-
-This command shows the host and device environments and software version. If PhasicFlow was build correctly, you would get the following output:
-```
-Initializing host/device execution spaces . . . 
-  Host execution space is Serial
-  Device execution space is Cuda
-
-  ou are using phasicFlow v-0.1 (copyright(C): www.cemf.ir)
-  In this build, double is used for floating point operations.
-
-
-Finalizing host/device execution space ....
-```

doc/mdDocs/phasicFlowFeatures.md

@@ -1,64 +1,116 @@
-# PhasicFlow Features {#phasicFlowFeatures}
+# PhasicFlow Features (v-1.0)
+
+The features of PhasicFlow described here are the main features that are implemented in the code for version 1.0. This document is not a complete list of all the features of PhasicFlow. The features are being added to the code continuously and this document may be behind the latest updates. Of course, the code review will give you the complete list.
+
+## Table of Contents
+
+- [1. Building options](#1-building-options)
+- [2. Preprocessing tools](#2-preprocessing-tools)
+- [3. Solvers for simulations](#3-solvers-for-simulations)
+- [4. Postprocessing tools](#4-postprocessing-tools)
+- [5. Models and features for simulations](#5-models-and-features-for-simulations)
+  - [5.1. General representation of walls](#51-general-representation-of-walls)
+  - [5.2. High precision integeration methods](#52-high-precision-integeration-methods)
+  - [5.3. Contact force models](#53-contact-force-models-needs-improvement)
+  - [5.4. Particle insertion](#54-particle-insertion)
+  - [5.5. Restarting/resuming a simulation](#55-restartingresuming-a-simulation)
+  - [5.6. Postprocessing data during simulation](#56-postprocessing-data-during-simulation)
+
+## 1. Building options
+
+You can build PhasicFlow to be executed on multi-core CPUs or GPUs. It is also possible to select the type of floating point variables in PhasicFlow: double or float. float type requires less memory and mostly consumes less time of a processor to complete a mathematical operation. So, there is a benefit for using floats in DEM simulation specially when GPU is targeted for computations.
 
-## Building options
-You can build PhasicFlow to be executed on multi-core CPUs or GPUs. It is also possible to select the type of floating point variables in PhasicFlow: double or float. float type requires less memory and mostly consumes less time of a processor to complete a mathematical operation. So, there is a benefit for using floats in DEM simulation specially when GPU is targeted for computations. 
 Build options for PhasicFlow:
-* **serial (double or float type)**: execution on one cpu core
-* **OpenMp (double or float type)**: execution on multiple cores of a CPU
-* **cuda (double or float type)**: execution on cuda-enabled GPUs
 
+- **serial (double or float type)**: execution on one cpu core
+- **OpenMp (double or float type)**: execution on multiple cores of a CPU
+- **cuda (double or float type)**: execution on cuda-enabled GPUs
 
+for more information on building PhasicFlow, please refer to the [installation guide](./howToBuild-V1.0.md).
 
-## Preprocessing tools
-Preprocessing tools are used to facilitate the process of case setup. They include tools for defining initial state of particles and geometry conversion. 
-* **particlesPhasicFlow** tool can be used to define the initial position of particles (for example at t = 0 s) and to set the initial field values for particles (like velocity, orientation, acceleration and etc).
-* **geometryPhasicFlow** converts user inputs for walls into a data structures that is used by PhasicFlow.
+## 2. Preprocessing tools
 
+PhasicFlow provides a set of tools for preprocessing the simulation case. These tools are used to define the initial state of particles, walls and other parameters that are required for running a simulation.
 
+- [**particlesPhasicFlow**](./../../utilities/particlesPhasicFlow/) tool can be used to define the initial position of particles (for example at t = 0 s) and to set the initial field values for particles (like velocity, orientation, acceleration, etc.).
 
-## Models and features for simulations
+- [**geometryPhasicFlow**](./../../utilities/geometryPhasicFlow/) converts user inputs for walls into a data structure that is used by PhasicFlow.
 
+## 3. Solvers for simulations
 
-### General representation of walls
+- [**sphereGranFlow**](./../../solvers/sphereGranFlow/) is a solver for simulating the flow of spherical particles with particle insertion mechanism. A full set of tutorial on various possible simulations can be found here: [sphereGranFlow tutorial](./../../tutorials/sphereGranFlow/).
+- [**grainGranFlow**](./../../solvers/grainGranFlow/) is a solver for simulating the flow of course-grained particles with particle insertion mechanism. A full set of tutorial on various possible simulations can be found here: [grainGranFlow tutorial](./../../tutorials/grainGranFlow/).
+- [**iterateGeometry**](./../../solvers/iterateGeometry/) is a  solver testing motion of walls without simulating particles. Since simulating with particles may take a long time and we may want to check the motion of geometry to be correct before actual simulation, we created this utility to test the motion of walls. A set of tutorial on various possible simulations can be found here: [iterateGeometry tutorial](./../../tutorials/iterateGeometry/).
+
+## 4. Postprocessing tools
+
+- [**pFlowToVTK**](./../../utilities/pFlowToVTK) is used to convert simulation results into vtk file format. vtk file format can be read by Paraview for visualizing the results.
+- [**postprocessPhasicFlow**](./../../utilities/postprocessPhasicFlow/) is a tool for performing various averaging and summation on the fields. Particle probing is also possible.
+
+## 5. Models and features for simulations
+
+### 5.1. General representation of walls
 
 Walls can be defined in three ways in PhasicFlow:
-* **Builtin walls** in PhasicFlow that include plane wall, cylinder/cone wall, cuboid, circle.
-* **stl wall**  that reads the data of the wall from an ASCII stl file.
-* **foamPatch wall** that reads the OpenFOAM mesh and converts the boundary patches into PhasicFlow walls (this feature is only available when performing CFD-DEM simulation using OpenFOAM).
 
-Walls can be fixed or in motion during simulations. Various motion models are implemented to cover most of the wall motions in phasicFlow ([see the source code] (./../../../src/MotionModel/)):
-* **fixedWall**  model, in which all walls are fixed. This model is mostly useful for granular flow under gravity or gas-solid flows (CFD-DEM). 
-* **rotatingAxisMotion**  model, in which walls are rotating around an axis of rotation with specified rotation speed. This model covers a wide range of granular flows in which the whole or a part of geometry is rotating, like mixers. 
-* **multiRotatingAxisMotion** model, in which a combination of rotations can be specified. One axis of rotation can itself have another axis of rotation, and so on. This creates the possibility of defining very complex motion pattern for walls, like what we see in Nauta blenders.
-* **vibratingMotion** model, in which walls vibrates based on a sinusoidal model with specified frequency and amplitude.
-In addition to these models, the user can add other motion models to the code based on their need. 
+- **Builtin walls** in PhasicFlow that include plane wall, cylinder/cone wall, cuboid, circle.
+- **stl wall**  that reads the data of the wall from an ASCII stl file.
+- **foamPatch wall** that reads the OpenFOAM mesh and converts the boundary patches into PhasicFlow walls (this feature is only available when performing CFD-DEM simulation using OpenFOAM).
 
+Walls can be fixed or in motion during simulations. Various motion models are implemented to cover most of the wall motions in phasicFlow ([see the source code](./../../src/MotionModel/)):
+
+- **stationay**  model, in which all walls are fixed. This model is mostly useful for granular flow under gravity or gas-solid flows (CFD-DEM).
+- **rotatingAxis**  model, in which walls are rotating around an axis of rotation with specified rotation speed. This model covers a wide range of granular flows in which the whole or a part of geometry is rotating, like mixers.
+- **multiRotatingAxis** model, in which a combination of rotations can be specified. One axis of rotation can itself have another axis of rotation, and so on. This creates the possibility of defining very complex motion pattern for walls, like what we see in Nauta blenders.
+- **vibrating** model, in which walls vibrates based on a sinusoidal model with specified frequency and amplitude.
+In addition to these models, the user can add other motion models to the code based on their need.
+
+### 5.2. High precision integeration methods
+
+The precision of integration in a DEM simulation is very important. Since sudden changes in the interaction forces occur during simulations (when objects contact or when they rebound). High precision integration methods makes it possible to accurately track position and velocity of objects (specially when they are in contact). When using these methods, it is possible to choose larger time steps for integration without loosing accuracy and causing instability in the simulation. Although  a high-precision integration requires more computations, but the benefits of choosing larger time steps in simulation can totally compensate it.
 
-### High precision integeration methods
-The precision of integration in a DEM simulation is very important. Since sudden changes in the interaction forces occur during simulations (when objects contact or when they rebound). High precision integration methods makes it possible to accurately track position and velocity of objects (specially when they are in contact). When using these methods, it is possible to choose larger time steps for integration without loosing accuracy and causing instability in the simulation. Although  a high-precision integration requires more computations, but the benefits of choosing larger time steps in simulation can totally compensate it. 
 Various integration methods are implemented in PhasicFlow:
 
-| Integration Method | Order | Type|
+|Integration Method | Order | Type|
 | :--- | :---: | :---: |
 | AdamsBashforth2 | 2 | one-step |
 | AdamsBashforth3 | 3 | one-step |
 | AdamsBashforth4 | 4 | one-step |
 | AdamsBashforth5 | 5 | one-step |
-| AdamsMoulton3 | 3 | predictor-corrector |
-| AdamsMoulton4 | 4 | predictor-corrector |
-| AdamsMoulton5 | 5 | predictor-corrector |
+| AdamsMoulton3 | 3 | predictor-corrector (not active)|
+| AdamsMoulton4 | 4 | predictor-corrector (not active)|
+| AdamsMoulton5 | 5 | predictor-corrector (not active)|
 
+### 5.3. Contact force models (needs improvement)
 
-### Contact force models
-Linear and non-linear visco-elastic contact force models are considered in the simulation. In addition to these, limited and non-limited Coulomb's friction model can be used to account for the friction between objects. For spherical objects, rolling friction can also be specified between bodies in contact. 
+Linear and non-linear visco-elastic contact force models are considered in the simulation. In addition to these, limited and non-limited Coulomb's friction model can be used to account for the friction between objects. For spherical objects, rolling friction can also be specified between bodies in contact.
+In addition, for course-grained particles simulation, we developed a speciall set of***
 
-### Particle insertion
-Particles can be inserted during simulation from specified region at specified rate and time interval. Any number of insertion regions can be defined in a simulation. Various region types are considered here: box, cylinder and sphere. Particles are inserted into the simulation through the specified region.
+### 5.4. Particle insertion
 
-### restarting/resuming a simulation
-It is possible to resume a simulation fron any time-folder that is avaiable in the simulation case setup directory. PhasicFlow restart the simulation from that time folder.
+Particles can be inserted during simulation from specified region at specified rate and time interval. Any number of insertion regions can be defined in a simulation. Various region types are considered here: `box`, `cylinder` and `sphere`. Particles are inserted into the simulation through the specified region.
 
-## Postprocessing tools
+### 5.5. restarting/resuming a simulation
 
-* **pFlowToVTK** is used to convert simulation results into vtk file format. vtk file format can be read by Paraview for visualizing the results.
-* **postprocessPhasicFlow** is a tool for performing various cell-based averaging on the fields. 
+It is possible to resume a simulation from any time-folder that is available in the simulation case setup directory. PhasicFlow restarts the simulation from that time folder.
+
+### 5.6. Postprocessing data during simulation
+
+PhasicFlow provides a powerful in-simulation postprocessing module that allows users to analyze particle data in real-time while the simulation is running. This feature enables:
+
+- **Real-time data analysis** without waiting for simulation completion
+- **Region-based processing** in spheres, along lines, or at specific points
+- **Various statistical operations** including weighted averages and sums of particle properties
+- **Individual particle tracking** to monitor specific particles throughout simulation
+- **Multiple processing methods** including arithmetic mean, uniform distribution, and Gaussian distribution
+- **Particle filtering** based on properties like diameter, mass, etc.
+- **Flexible time control** options for when postprocessing should be executed
+
+To activate in-simulation postprocessing, users need to:
+
+1. Create a `postprocessDataDict` file in the `settings` directory with appropriate configurations
+2. Add `libs ("libPostprocessData.so")` and `auxFunctions postprocessData` to the `settings/settingsDict` file
+
+Results are written to output files in the case directory with timestamps, allowing users to monitor simulation behavior as it progresses without interrupting the simulation. for more information on how to use this feature, please refer to the [PostprocessData](./../../src/PostprocessData/) module.
+
+The same postprocessing module can also be used after simulation completion through the [`postprocessPhasicFlow`](./../../utilities/postprocessPhasicFlow/) utility.

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

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

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

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

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

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

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

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

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

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

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

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

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

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

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

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

src/Integration/integration/integration.cpp

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

src/Integration/integration/integration.hpp

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

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

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

src/Particles/GrainParticles/grainShape/grainShape.cpp

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

src/Particles/GrainParticles/grainShape/grainShape.hpp

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

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

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

src/Particles/SphereParticles/sphereShape/sphereShape.cpp

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

src/Particles/SphereParticles/sphereShape/sphereShape.hpp

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

src/Particles/dynamicPointStructure/dynamicPointStructure.cpp

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

src/Particles/particles/particles.hpp

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

src/Particles/particles/shape/shape.cpp

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

src/Particles/particles/shape/shape.hpp

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

src/PostprocessData/CMakeLists.txt

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

src/PostprocessData/fieldsDataBase/fieldFunctions.hpp

@@ -24,7 +24,7 @@ Licence:
 #include "types.hpp"
 #include "span.hpp"
 
-namespace pFlow
+namespace pFlow::postprocessData
 {
 
 template<typename T>

src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp

@@ -21,31 +21,45 @@ 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)
+
+bool pFlow::postprocessData::fieldsDataBase::loadPointStructureToTime()
+{
+    return false;
+}
+
+
+bool pFlow::postprocessData::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;
 }
 
-pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetRealField(const word &name)
+pFlow::span<pFlow::real> pFlow::postprocessData::fieldsDataBase::createOrGetRealField
+(
+    const word &name
+)
 {
     
     bool shouldUpdate = checkForUpdate(name);
@@ -70,7 +84,181 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetRealField(const word
         field.size());
 }
 
-bool pFlow::fieldsDataBase::findFunction(
+pFlow::span<pFlow::real> pFlow::postprocessData::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::postprocessData::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::postprocessData::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::postprocessData::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::postprocessData::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::postprocessData::fieldsDataBase::findFunction(
     const word &compoundFieldName,
     word &fieldName,
     fieldsDataBase::Functions &func)
@@ -176,7 +364,7 @@ bool pFlow::fieldsDataBase::findFunction(
     return false; // No match
 }
 
-bool pFlow::fieldsDataBase::inputOutputType
+bool pFlow::postprocessData::fieldsDataBase::inputOutputType
 (
     fieldsDataBase::Functions func, 
     const word &inputType, 
@@ -274,95 +462,145 @@ bool pFlow::fieldsDataBase::inputOutputType
     return false;
 }
 
-pFlow::fieldsDataBase::fieldsDataBase(const Time &time, bool inSimulation)
-: 
-    time_(time), 
+pFlow::postprocessData::fieldsDataBase::fieldsDataBase
+(
+    systemControl& control,
+    const dictionary& postDict,
+    bool inSimulation,
+    timeValue startTime
+)
+:
+    time_(control.time()),
     inSimulation_(inSimulation)
-{}
+{
+    if(!inSimulation_)
+    {
+        shapeType_ = postDict.getValOrSet<word>("shapeType", "");
 
-pFlow::timeValue pFlow::fieldsDataBase::currentTime() const
+        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::postprocessData::fieldsDataBase::currentTime() const
 {
     return time_.currentTime();
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::postprocessData::fieldsDataBase::getFieldTypeNameFunction
 (
-    const word& compoundName, 
-    word& fieldName,
+    const word& compoundName,
+    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::postprocessData::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::postprocessData::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)
+pFlow::span<pFlow::realx3> pFlow::postprocessData::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(), 
@@ -370,7 +608,7 @@ pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updatePoints(bool forceUpdate)
 
 }
 
-pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updateFieldRealx3
+pFlow::span<pFlow::realx3> pFlow::postprocessData::fieldsDataBase::updateFieldRealx3
 (
     const word &compoundName, 
     bool forceUpdate
@@ -379,7 +617,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;
@@ -401,7 +644,7 @@ pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updateFieldRealx3
         
 }
 
-pFlow::span<pFlow::realx4> pFlow::fieldsDataBase::updateFieldRealx4
+pFlow::span<pFlow::realx4> pFlow::postprocessData::fieldsDataBase::updateFieldRealx4
 (
     const word &compoundName, 
     bool forceUpdate
@@ -410,7 +653,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;
@@ -432,7 +680,7 @@ pFlow::span<pFlow::realx4> pFlow::fieldsDataBase::updateFieldRealx4
     
 }
 
-pFlow::span<pFlow::real> pFlow::fieldsDataBase::updateFieldReal
+pFlow::span<pFlow::real> pFlow::postprocessData::fieldsDataBase::updateFieldReal
 (
     const word &compoundName, 
     bool forceUpdate
@@ -441,7 +689,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;
@@ -594,7 +847,7 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::updateFieldReal
     return span<real>(nullptr, 0);
 }
 
-pFlow::span<pFlow::uint32> pFlow::fieldsDataBase::updateFieldUint32
+pFlow::span<pFlow::uint32> pFlow::postprocessData::fieldsDataBase::updateFieldUint32
 (
     const word& name,
     bool forceUpdate
@@ -603,17 +856,15 @@ pFlow::span<pFlow::uint32> pFlow::fieldsDataBase::updateFieldUint32
     return updateField<uint32>(name, forceUpdate);
 }
 
-pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
+pFlow::postprocessData::allPointFieldTypes pFlow::postprocessData::fieldsDataBase::updateFieldAll
 (
     const word &compoundName, 
     bool forceUpdate
 )
 {
-
-    word originalType, typeAfterFunction, fieldName;
-    Functions func;
+    word typeAfterFunction, originalType;
     
-    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 +872,7 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
         return span<real>(nullptr, 0);
     }
 
-
-    if(  typeAfterFunction== getTypeName<realx3>()  )
+    if(  typeAfterFunction == getTypeName<realx3>()  )
     {
         return updateFieldRealx3(compoundName, forceUpdate);
     }
@@ -642,22 +892,12 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
     }
 }
 
-const pFlow::pointStructure &pFlow::fieldsDataBase::pStruct() const
-{
-    if(inSimulation_)
-    {
-        return 
-            static_cast<const pointStructure&>(
-                time_.lookupObject<dynamicPointStructure>(pointStructureFile__)
-            );
-    }
-    else
-    {
-        return time_.lookupObject<pointStructure>(pointStructureFile__);
-    }
-}
-
-bool pFlow::pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace)
+bool pFlow::postprocessData::fieldsDataBase::pointFieldGetType
+(
+    const word& TYPENAME, 
+    word& fieldType, 
+    word& fieldSpace
+)
 {
 	std::regex match("pointField\\<([A-Za-z1-9_]*)\\,([A-Za-z1-9_]*)\\>");
     std::smatch search;
@@ -668,3 +908,46 @@ bool pFlow::pointFieldGetType(const word& TYPENAME, word& fieldType, word& field
 	return true;
 }
 
+pFlow::uniquePtr<pFlow::postprocessData::fieldsDataBase> 
+    pFlow::postprocessData::fieldsDataBase::create
+(
+    systemControl& control, 
+    const dictionary& postDict, 
+    bool inSimulation,
+    timeValue startTime
+)
+{
+    word dbType;
+    if(inSimulation)
+    {
+        dbType = "fieldsDataBase<simulation>";
+    }
+    else
+    {
+        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;
+}
+
+
+

src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp

@@ -21,52 +21,26 @@ 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;
+}
 
-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>>; 
+namespace pFlow::postprocessData
+{
 
 
 class fieldsDataBase
-
 {
 private:
 
@@ -99,9 +73,6 @@ private:
 
 
     // - Member variables
-    
-        /// const reference to the Time object
-        const Time&         time_;
 
         /// List to store all the point fields 
         anyList             allFields_;
@@ -109,114 +80,105 @@ private:
         /// Map to store the last capture time of each field
         wordMap<timeValue>  captureTime_;
 
-        /// 
-        bool                inSimulation_ = false;
-        
-        static const inline std::map<word, word> reservedFieldNames_
-        {
-            {"position", "realx3"},
-            {"one", "real"}
-        };
-        
-    // - Private member functions
+        /// Reference to the Time object
+        Time&               time_; 
 
+        /// Flag indicating if we're in simulation mode
+        bool                inSimulation_ = false;
+
+        word                shapeType_;
+
+protected:
+
+        /// check if pointField name exists in Time or time folder 
+        virtual 
+        bool pointFieldNameExists(const word& name)const = 0;
+        
+        /// Loads a pointField with a given name to the Time object 
+        virtual 
+        bool loadPointFieldToTime(const word& name)= 0;
+
+        virtual 
+        bool loadPointStructureToTime()=0;
+        
+        virtual 
+        const shape& getShape() const= 0;
+
+        const word& shapeTypeName()const
+        {
+            return shapeType_;
+        }
+
+        /// get the type name of the pointField in the Time object 
+        virtual 
+        word getPointFieldType(const word& name)const = 0;
+
+        /// Checks if a field needs to be updated based on capture time
+        bool checkForUpdate(const word& compoundName, bool forceUpdate = false);
+    
+        /// @brief return the size of pointStructure 
         uint32 pointFieldSize()
         {
             auto s = updatePoints();
             return s.size();
         }
-        
-        /// @brief Checks if a field needs to be updated.
-        /// @param name The name of the field to check.
-        /// @param forceUpdate Forces an update if true. Defaults to `false`.
-        /// @return `true` if the field needs updating or `forceUpdate` is true, `false` otherwise.
-        bool checkForUpdate(const word& name, bool forceUpdate = false);
 
-        /**
-         * @brief Updates and retrieves a point field of a specified type from the database.
-         * 
-         * This is a template function that updates and retrieves a point field of type `T`
-         * from the database. It checks if the field needs to be updated based on the last
-         * capture time or if a forced update is requested. If an update is necessary, it
-         * retrieves the latest data for the field.
-         * 
-         * @tparam T The type of the point field to update and retrieve. Must be a ValidFieldType.
-         * @param name The name of the field to update.
-         * @param forceUpdate A boolean flag indicating whether to force an update of the field
-         *                    regardless of its current state. Defaults to `false`.
-         * 
-         * @return A span of `T` representing the updated field data.
-         * 
-         * @throws message If the field cannot be found in the database or if there is a type mismatch.
-         */
         template<ValidFieldType T>
         span<T> updateField(const word& name, bool forceUpdate = false);
+
+        template<ValidFieldType T>
+        span<T> updateReservedField(const word& name, bool forceUpdate = false);
         
-        /// @brief Creates a new real field or retrieves an existing one.
-        /// 
-        /// If a field with the given name already exists, it returns a span to that field.
-        /// If the field does not exist, it creates a new real field with the given name
-        /// and returns a span to the newly created field.
-        /// 
-        /// @param name The name of the field to create or retrieve.
-        /// @return span<real> of the field 
         span<real> createOrGetRealField(const word& name);
 
-        /**
-         * @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,140 +193,75 @@ public:
         fieldsDataBase& operator=(fieldsDataBase&&) = delete;
         
         /// destructor
-        ~fieldsDataBase() = default;
+        virtual ~fieldsDataBase() = default;
+
+        create_vCtor
+        (
+            fieldsDataBase,
+            bool,
+            (
+                systemControl& control, 
+                const dictionary& postDict, 
+                bool inSimulation,
+                timeValue startTime
+            ),
+            (control, postDict, inSimulation, startTime)
+        );
+  
     
     // - Public Access Functions
         /// returns the current time
         timeValue currentTime()const;
 
-        /// const ref to object Time
-        const Time& time()const
+        /// const ref to object Time 
+        const Time& time()const 
         {
             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);   
-        
-        /// overload for the function getPointFieldType without `func` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& originalType, 
-            word& typeAfterFunction);
-        
-        /// overload for function getPointFieldType without `originalType` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& typeAfterFunction);
+            Functions& func
+        )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.
-         */
+        
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& originalType, 
+            word& typeAfterFunction
+        ) const;
+
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& typeAfterFunction
+        ) const;
+
+        
+        /// 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);
@@ -378,11 +275,51 @@ public:
         allPointFieldTypes updateFieldAll(
             const word& compoundName, 
             bool forceUpdate = false);
+        
+        virtual 
+        const pointStructure& pStruct()const = 0;
+        
+        /// Get the next avaiable time folder after the current time folder 
+        /// This is only used for post-simulation processing 
+        virtual
+        timeValue getNextTimeFolder()const
+        {
+            return -1.0;
+        }
 
-        const pointStructure& pStruct()const;
+        /// Sets the current folder to the next time folder.
+        /// This is used only for post-simulation processing
+        /// @returns the time value of the next folder.
+        virtual 
+        timeValue setToNextTimeFolder()
+        {
+            return -1.0;
+        }
+        
+        /// Skips the next time folder.
+        /// This is used only for post-simulation processing 
+        /// @returns the time value of the skipped folder 
+        virtual 
+        timeValue skipNextTimeFolder()
+        {
+            return -1.0;
+        }
+    
+    static
+    bool pointFieldGetType(
+        const word& TYPENAME, 
+        word& fieldType, 
+        word& fieldSpace);
+
+    static
+    uniquePtr<fieldsDataBase> create(
+        systemControl& control, 
+        const dictionary& postDict, 
+        bool inSimulation,
+        timeValue startTime);
 };
 
-} // nampespace pFlow
+} // namespace pFlow::postprocessData
 
 #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::postprocessData
+{
+
+
+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

@@ -23,49 +23,34 @@ Licence:
 
 #include "fieldsDataBase.hpp"
 
-template<pFlow::ValidFieldType T>
+template<pFlow::postprocessData::ValidFieldType T>
 inline
-pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUpdate)
+pFlow::span<T> pFlow::postprocessData::fieldsDataBase::updateField(const word& name, bool forceUpdate)
 {
   
     bool shouldUpdate = checkForUpdate(name, forceUpdate);
         
     if(shouldUpdate)
     {
-        if(name == "one")
-        {   
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>
-            (
-                "one",
-                FieldTypeHost<T>
-                (
-                    "one", 
-                    "value",
-                    pointFieldSize(),
-                    T{1}
-                )
-            );
-        }
-        else if( name == "position")
+        if(reservedFieldNames_.contains(name))
         {
-            if constexpr( std::same_as<T, realx3>)
+            return updateReservedField<T>(name, true);
+        }   
+        else
+        {
+            if( loadPointFieldToTime(name) )
             {
-                return updatePoints(forceUpdate);            
+                const auto& pField = time_.template lookupObject<pointField_D<T>>(name);
+                allFields_.template emplaceBackOrReplace<FieldTypeHost<T>>(
+                name, 
+                pField.activeValuesHost());
             }
             else
             {
-                fatalErrorInFunction<< "Error in getting the type name of field: "<<
-                    name<<", with type name: "<< getTypeName<T>() <<endl;
+                fatalErrorInFunction<<"Error in loading the pointField "<<name<<endl;
                 fatalExit;
-                return span<T>(nullptr, 0);
             }
-        }
-        else
-        {
-            const auto& pField = time_.lookupObject<pointField_D<T>>(name);
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
-            name, 
-            pField.activeValuesHost());
+            
         }
     }
 
@@ -77,4 +62,110 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
     
 }
 
+
+template<pFlow::postprocessData::ValidFieldType T>
+inline
+pFlow::span<T> pFlow::postprocessData::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"
+
+
+
+bool pFlow::postprocessData::simulationFieldsDataBase::pointFieldNameExists(const word &name) const
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::postprocessData::simulationFieldsDataBase::loadPointFieldToTime(const word &name)
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::postprocessData::simulationFieldsDataBase::loadPointStructureToTime()
+{
+    // it is already in the Time object 
+    return time().lookupObjectName(pointStructureFile__);
+}
+
+
+const pFlow::shape& pFlow::postprocessData::simulationFieldsDataBase::getShape() const
+{
+    return shape_;
+}
+
+pFlow::word pFlow::postprocessData::simulationFieldsDataBase::getPointFieldType
+(
+    const word &name
+) const
+{
+    word pfType =  time().lookupObjectTypeName(name);
+    word type, space;
+    if(!fieldsDataBase::pointFieldGetType(pfType, type, space))
+    {
+        fatalErrorInFunction
+            <<"Error in retriving the type of pointField "
+            << pfType<<endl;
+        fatalExit;
+    }
+    return type;
+}
+
+pFlow::postprocessData::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::postprocessData::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::postprocessData
+{
+
+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,25 @@
+#include "PostprocessOperationAvMassVelocity.hpp"
+
+namespace pFlow::postprocessData
+{
+
+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::postprocessData
+{
+
+
+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
+    );
+
+}; 
+
+
+} // namespace pFlow::postprocessData
+
+#endif //__PostprocessOperationAvMassVelocity_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -0,0 +1,181 @@
+#include "PostprocessOperationAverage.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+#include "operationFunctions.hpp"
+
+namespace pFlow::postprocessData
+{
+
+/// Constructs average processor and initializes result field based on input field type
+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;
+    }
+}
+
+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 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;
+}
+
+bool PostprocessOperationAverage::write(const fileSystem &parDir) const
+{   
+    if(! postprocessOperation::write(parDir))
+    {
+        return false;
+    }
+    if(!calculateFluctuation2_())
+    {
+        return true;
+    }
+    
+    auto ti = time().TimeInfo();
+
+    if(!os2Ptr_)
+    {
+        fileSystem path = parDir+(
+            processedFieldName()+"_prime2" + ".Start_" + ti.timeName());
+        os2Ptr_ = makeUnique<oFstream>(path);
+        
+        regPoints().write(os2Ptr_());
+    }
+
+    
+    std::visit
+    (
+        [&](auto&& arg)->bool
+        {
+            return writeField(os2Ptr_(), ti.t(), arg, threshold());
+        },
+        fluctuation2FieldPtr_()
+    );
+
+    return true;
+}
+
+} // namespace pFlow::postprocessData

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

@@ -0,0 +1,211 @@
+/*------------------------------- 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::postprocessData
+{
+
+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;
+
+    /// Pointer to the output stream for writing fluctuation2 results
+    mutable uniquePtr<oFstream>       os2Ptr_ = 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_();
+    }
+
+    /// write to os stream
+    bool write(const fileSystem &parDir)const override;
+
+
+    /// @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;
+
+}; 
+
+
+} // namespace pFlow::postprocessData
+
+#endif //__PostprocessOperationAverage_hpp__