Merge pull request #198 from PhasicFlow/postProcessing

Post processing

src/PostprocessData/CMakeLists.txt

@@ -10,6 +10,7 @@ set(SourceFiles
     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 +20,7 @@ set(SourceFiles
     # Operations
     operation/postprocessOperation/postprocessOperation.cpp
     operation/PostprocessOperation/PostprocessOperationSum.cpp
+    operation/PostprocessOperation/PostprocessOperationAverage.cpp
     operation/includeMask/includeMask.cpp
     operation/includeMask/IncludeMasks.cpp
 

src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp

@@ -33,7 +33,7 @@ bool pFlow::fieldsDataBase::checkForUpdate(const word &name, bool forceUpdate)
     
     if(auto [iter, found]= captureTime_.findIf(name); found)
     {
-        shouldUpdate = iter->second < t;
+        shouldUpdate = iter->second < t || forceUpdate;
         iter->second = t;
     }
     else

src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp

@@ -50,7 +50,7 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
         {
             if constexpr( std::same_as<T, realx3>)
             {
-                return updatePoints(forceUpdate);            
+                return updatePoints(true);            
             }
             else
             {

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp

@@ -0,0 +1,173 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __PostprocessOperationAvMassVelocity_hpp__
+#define __PostprocessOperationAvMassVelocity_hpp__
+
+/*!
+ * @class PostprocessOperationAvMassVelocity
+ * @brief A class for averaging field values within specified regions during post-processing.
+ * 
+ * @details
+ * The PostprocessOperationAvMassVelocity class is a specialized post-processing operation that 
+ * calculates the average of field values within specified regions. It inherits from the 
+ * postprocessOperation base class and implements a weighted averaging operation that 
+ * can be applied to scalar (real), vector (realx3), and tensor (realx4) fields.
+ * 
+ * The average operation follows the mathematical formula:
+ * \f[
+ *    \text{result} = \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * Where:
+ * - \f$ i \f$ represents all particles within the specified processing region
+ * - \f$ j \f$ belongs to a subset of \f$ i \f$ based on an includeMask
+ * - \f$ w_i \f$ is the weight factor for particle \f$ i \f$
+ * - \f$ \phi_i \f$ is the value from the phi field for particle \f$ i \f$
+ * - \f$ \text{field}_j \f$ is the value from the target field for particle \f$ j \f$
+ * 
+ * The calculation can optionally be divided by the region volume (when divideByVolume is set to yes),
+ * which allows calculating normalized averages:
+ * \f[
+ *    \text{result} = \frac{1}{V_{\text{region}}} \frac{\sum_{j \in \text{includeMask}} w_j \cdot \phi_j \cdot \text{field}_j}
+ *                        {\sum_{i \in \text{processRegion}} w_i \cdot \phi_i}
+ * \f]
+ * 
+ * The averaging can be further filtered using an includeMask to selectively include only
+ * specific particles that satisfy certain criteria.
+ * 
+ * This class supports the following field types:
+ * - real (scalar values)
+ * - realx3 (vector values)
+ * - realx4 (tensor values)
+ * 
+ * @section usage Usage Example
+ * Below is a sample dictionary showing how to configure and use this class:
+ *
+ * ```
+ * processMethod    arithmetic;     // method of performing the sum (arithmetic, uniformDistribution, GaussianDistribution)
+ * processRegion    sphere;         // type of region on which processing is performed 
+ * 
+ * sphereInfo
+ * {
+ *     radius 0.01;
+ *     center (-0.08 -0.08 0.015);
+ * }
+ * 
+ * timeControl default; 
+ * 
+ * /// all the post process operations to be done 
+ * operations
+ * (
+ *     // computes the arithmetic mean of particle velocity 
+ *     averageVel
+ *     {
+ *         function         average;
+ *         field            velocity;
+ *         dividedByVolume  no;     // default is no
+ *         threshold        3;      // default is 1
+ *         includeMask      all;    // include all particles in the calculation
+ *     }
+ *     
+ *     // computes the fraction of par1 in the region 
+ *     par1Fraction
+ *     {
+ *         function        average;
+ *         field           one;     // the "one" field is special - all members have value 1.0
+ *         phi             one;     // default is "one"
+ *         dividedByVolume no;
+ *         includeMask     lessThan;
+ *         
+ *         // diameter of par1 is 0.003, so these settings 
+ *         // will select only particles of type par1
+ *         lessThanInfo
+ *         {
+ *             field      diameter;
+ *             value      0.0031;  
+ *         }
+ *     }
+ * );
+ * ```
+ *
+ * @section defaults Default Behavior
+ * - By default, `phi` is set to the field named "one" which contains value 1.0 for all entries
+ * - `dividedByVolume` is set to "no" by default
+ * - `threshold` is set to 1 by default
+ * - `includeMask` can be set to various filters, with "all" being the default to include all particles
+ *
+ * @section special Special Fields
+ * The field named "one" is a special field where all members have the value 1.0. This makes it
+ * particularly useful for calculating:
+ * 
+ * 1. Volume or number fractions (as shown in the par1Fraction example)
+ * 2. Simple counts when used with an appropriate mask
+ * 3. Normalizing values by particle count
+ *
+ * @see postprocessOperation
+ * @see executeAverageOperation
+ */
+
+#include <variant>
+#include <vector>
+
+#include "postprocessOperation.hpp"
+#include "regionField.hpp"
+#include "includeMask.hpp"
+
+namespace pFlow
+{
+
+
+class PostprocessOperationAvMassVelocity
+:
+    public postprocessOperation
+{
+   
+public:
+
+    TypeInfo("PostprocessOperation<avMassVelocity>");
+
+    /// @brief Constructs average operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
+    PostprocessOperationAvMassVelocity(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    /// destructor 
+    ~PostprocessOperationAvMassVelocity() override = default;
+
+    /// add this virtual constructor to the base class
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationAvMassVelocity,
+        dictionary
+    );
+
+}; 
+
+
+}
+
+#endif //__PostprocessOperationAvMassVelocity_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -0,0 +1,138 @@
+#include "PostprocessOperationAverage.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+#include "fieldFunctions.hpp"
+
+/// Constructs average processor and initializes result field based on input field type
+pFlow::PostprocessOperationAverage::PostprocessOperationAverage
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(opDict, regPoints, fieldsDB),
+    calculateFluctuation2_(opDict.getValOrSet<Logical>("fluctuation2", Logical(false)))
+{
+    if( fieldType() == getTypeName<real>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, real(0)));
+    }
+    else if( fieldType() == getTypeName<realx3>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx3(0)));
+    }
+    else if( fieldType() == getTypeName<realx4>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx4(0)));
+    }
+    else
+    {
+        fatalErrorInFunction<<" in dictionary "<< opDict.globalName()
+            << " field type is not supported for average operation"
+            << " field type is "<< fieldType()
+            << endl;
+        fatalExit;
+    }
+}
+
+pFlow::PostprocessOperationAverage::PostprocessOperationAverage
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName,
+    const word &includeName, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(opDict, fieldName, phiName, includeName, regPoints, fieldsDB),
+    calculateFluctuation2_(opDict.getValOrSet<Logical>("fluctuation2", Logical(false)))
+{
+    if( fieldType() == getTypeName<real>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, real(0)));
+    }
+    else if( fieldType() == getTypeName<realx3>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx3(0)));
+    }
+    else if( fieldType() == getTypeName<realx4>() )
+    {
+        processedRegFieldPtr_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx4(0)));
+    }
+    else
+    {
+        fatalErrorInFunction<<" in dictionary "<< opDict.globalName()
+            << " field type is not supported for average operation"
+            << " field type is "<< fieldType()
+            << endl;
+        fatalExit;
+    }
+}
+
+/// Performs weighted average of field values within each region
+bool pFlow::PostprocessOperationAverage::execute
+(
+    const std::vector<span<real>>& weights
+)
+{
+    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,
+                regP,
+                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),
+                        dbVol,
+                        weights,
+                        mask);
+                }
+            },
+            allField)
+        );
+    }
+    
+
+    return true;
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp

@@ -3,6 +3,7 @@
 #include "fieldsDataBase.hpp"
 #include "fieldFunctions.hpp"
 
+/// Constructs sum processor and initializes result field based on input field type
 pFlow::PostprocessOperationSum::PostprocessOperationSum
 (
     const dictionary &opDict, 
@@ -37,6 +38,7 @@ pFlow::PostprocessOperationSum::PostprocessOperationSum
     }
 }
 
+/// Performs weighted sum of field values within each region
 bool pFlow::PostprocessOperationSum::execute
 (
     const std::vector<span<real>>& weights

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp

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

src/PostprocessData/operation/PostprocessOperation/fieldFunctions.hpp

@@ -31,7 +31,6 @@ Licence:
 namespace pFlow
 {
 
-
 template<typename T>
 regionField<T> executeSumOperation
 (
@@ -45,13 +44,14 @@ regionField<T> executeSumOperation
 )
 {
     regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)
     {
         auto partIndices = regPoints.indices(reg);
-        auto vols = regPoints.volumes();
+        
         auto w = weights[reg];
-        T sum{};
+        T sum = T{};
         uint n = 0;
         for(auto index:partIndices)
         {
@@ -80,33 +80,106 @@ regionField<T> executeAverageOperation
 (
     const word&                    regFieldName,
     const span<T>&                 field, 
-    const regionPoints&            regPoints, 
+    const regionPoints&            regPoints,
+    const bool                     devideByVol, 
     const std::vector<span<real>>& weights,
     const span<real>&              phi,
     const includeMask::Mask&       mask
 )
 {
     regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
 
     for(uint32 reg =0; reg<regPoints.size(); reg++)
     {
         auto partIndices = regPoints.indices(reg);
         auto w = weights[reg];
-        T sumNum{};
+        T sumNum = T{};
-        real sumDen{};
+        real sumDen = 0;
         uint n = 0;
         for(auto index:partIndices)
         {
-            if( index!= -1 && mask( index ))
+            if( index!= -1)
             {
-                sumNum += w[n] * field[index]* phi[index];
+                if( mask(index))
+                {
+                    sumNum += w[n] * field[index]* phi[index];
+                }
+                sumDen += w[n] * phi[index];
             }
-            sumDen += w[n] * phi[index];
+            
             n++;
         }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / vols[reg];
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
+        
+    }
 
-        sumDen = max(sumDen, smallValue);
+    return processedField;
-        processedField[reg] = sumNum/sumDen;
+}
+
+
+template<typename T>
+regionField<T> executeFluctuation2Operation
+(
+    const word&                    regFieldName,
+    const span<T>&                 field,
+    const regionField<T>&          fieldAvg, 
+    const bool                     devideByVol, 
+    const std::vector<span<real>>& weights,
+    const includeMask::Mask&       mask
+)
+{
+    const auto& regPoints = fieldAvg.regPoints();
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto w = weights[reg];
+        auto vol = vols[reg];
+        T avField{};
+        if(devideByVol)
+        {
+            avField = vol * fieldAvg[reg];
+        }
+        else
+        {
+            avField = fieldAvg[reg];
+        }
+        
+        T sumNum = T{};
+        real sumDen = 0;
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( index!= -1)
+            {
+                if( mask(index))
+                {
+                    sumNum += w[n] * pow( avField- field[index],static_cast<real>(2));
+                }
+                sumDen += w[n];
+            }
+            n++;
+        }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / vol;
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
         
     }
 

src/PostprocessData/operation/includeMask/includeMask.cpp

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

src/PostprocessData/operation/includeMask/includeMask.hpp

@@ -74,6 +74,8 @@ public:
 
 	includeMask(const dictionary& opDict, fieldsDataBase& feildsDB);
 
+    includeMask(const word& type, const dictionary& opDict, fieldsDataBase& feildsDB);
+
 	virtual ~includeMask() = default;
 
     create_vCtor
@@ -85,6 +87,18 @@ public:
         ),
         (opDict, feildsDB)		
     );
+
+    create_vCtor
+    (
+        includeMask, 
+        word,
+        (
+            const word& type, 
+            const dictionary& opDict, 
+            fieldsDataBase& feildsDB
+        ),
+        (type, opDict, feildsDB)
+    );
 	 
 	const fieldsDataBase& database()const
 	{
@@ -103,6 +117,12 @@ public:
 	uniquePtr<includeMask> create(
         const dictionary& opDict, 
         fieldsDataBase& feildsDB);
+    
+    static
+    uniquePtr<includeMask> create(
+        const word& type,
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB);
 	
 };
 

src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp

@@ -82,14 +82,33 @@ bool writeField
 
 }
 
-
-
 pFlow::postprocessOperation::postprocessOperation
 (
     const dictionary &opDict, 
     const regionPoints& regPoints,
     fieldsDataBase &fieldsDB
 )
+:
+    postprocessOperation
+    (
+        opDict,
+        opDict.getVal<word>("field"),
+        opDict.getValOrSet<word>("phi", "one"),
+        opDict.getValOrSet<word>("includeMask", "all"),
+        regPoints,
+        fieldsDB
+    )
+{}
+
+pFlow::postprocessOperation::postprocessOperation
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName,
+    const word& includeName,  
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
 :
     operationDict_(opDict),
     threshold_
@@ -110,15 +129,15 @@ pFlow::postprocessOperation::postprocessOperation
     ),
     fieldName_
     (
-        opDict.getValOrSet<word>("field", "one")
+        fieldName
     ),
     phiFieldName_
     (
-        opDict.getValOrSet<word>("phi", "one")
+        phiName
     ),
     includeMask_
     (
-        includeMask::create(opDict, fieldsDB)
+        includeMask::create(includeName, opDict, fieldsDB)
     )
 {
 
@@ -128,7 +147,6 @@ pFlow::postprocessOperation::postprocessOperation
         fatalExit;
     }
 }
-
 const pFlow::Time& pFlow::postprocessOperation::time() const
 {
     return database_.time();

src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp

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

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -67,6 +67,10 @@ bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::execute
         return true;
     }
 
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
     // update processing methods
     auto& regPoints = this->regPoints();
     

src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp

@@ -92,7 +92,12 @@ bool pFlow::particleProbePostprocessComponent::execute
         executed_ = false;
         return true;
     }
-    
+
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
+
     if(!regionPointsPtr_().update())
     {
         fatalErrorInFunction

src/PostprocessData/postprocessData/postprocessData.cpp

@@ -40,14 +40,16 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
             objectFile::READ_IF_PRESENT,
             objectFile::WRITE_NEVER
         )
-    ),
+    )
-    componentsDicts_(readDictList("components", dict_))
 {
     postProcessGlobals::defaultDir__ = CWD()/pFlow::postProcessGlobals::defaultRelDir__;
     
     // if dictionary is not provided, no extra action is required. 
-    if( !dict_.fileExist() )
+    if( !dict_.fileExist() || !dict_.headerOk() )
     {
+        WARNING<<"You requested postprocessData function while,"
+            <<" the dictionary system/postprocessDataDict does not exist."
+            <<" This feature is disabled in the current run."<<END_WARNING;
         return;
     }
 
@@ -72,7 +74,9 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
             "execution");
     }
     
-    for(auto& compDict:componentsDicts_)
+    componentsDictsPtr_ = makeUnique<dictionaryList>(readDictList("components", dict_));
+
+    for(auto& compDict:*componentsDictsPtr_)
     {
         postprocesses_.push_back( postprocessComponent::create(
             compDict, 

src/PostprocessData/postprocessData/postprocessData.hpp

@@ -63,7 +63,7 @@ class postprocessData
     fileDictionary        dict_;
 
     /// list of dictionaries for postprocess components 
-    dictionaryList        componentsDicts_;
+    uniquePtr<dictionaryList>        componentsDictsPtr_ = nullptr;
 
     /// @brief  default time control that can be used for all post-process components
     uniquePtr<baseTimeControl> defaultTimeControlPtr_= nullptr;

src/PostprocessData/region/regionFields/regionField.hpp

@@ -31,6 +31,22 @@ namespace pFlow
 template<typename T> 
 class regionField
 {
+public:
+    using FieldType         = Field<T, HostSpace>;
+
+    using iterator        	= typename FieldType::iterator;
+
+    using const_iterator   	= typename FieldType::const_iterator;
+
+    using reference       	= typename FieldType::reference;
+    
+    using const_reference  	= typename FieldType::const_reference;
+
+    using value_type       	= typename FieldType::value_type;
+    
+    using pointer         	= typename FieldType::pointer;
+    
+    using const_pointer    	= typename FieldType::const_pointer;
 private:
            
     /// the field value

src/PostprocessData/region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.hpp

@@ -26,6 +26,34 @@ Licence:
 namespace pFlow
 {
 
+/**
+ * @class centerPointsRegionPoints
+ * @brief A region points implementation that selects particles based on their IDs
+ * 
+ * This class is responsible for selecting points (particles) by their IDs from
+ * a simulation database and tracking their properties. It maintains information 
+ * about the selected particles including their positions, volumes, and diameters.
+ * 
+ * The selection is performed based on IDs provided in the input dictionary.
+ * Once selected, the particles' properties can be accessed through various 
+ * methods. The update method allows refreshing the selection when particle data
+ * changes. The selection occurs at startTime defined in the time control, and
+ * there are some methods for selecting ids:
+ *  - specifying ids
+ *  - using selectors specified in pStructSelector class, which includes:
+ *    - box: selects particles within a box region
+ *    - sphere: selects particles within a spherical region
+ *    - cylinder: selects particles within a cylindrical region
+ *    - random: randomly selects a specified number of particles
+ *    - strided: selects particles with a specified stride pattern
+ * 
+ * This class is useful for tracking specific particles of interest throughout
+ * a simulation and analyzing their behavior.
+ * 
+ * @see regionPoints Base class providing the interface for different region 
+ *                   point selections
+ * @see pStructSelector Class providing different particle selection methods
+ */
 class centerPointsRegionPoints
 :
     public regionPoints
@@ -59,6 +87,7 @@ private:
 
 public:
 
+    /// Type info
     TypeInfo("centerPoints");
 
     centerPointsRegionPoints(
@@ -67,50 +96,69 @@ public:
 
     ~centerPointsRegionPoints() override = default;
 
+    /// @brief  Returns the number of selected points/particles
+    /// @return Number of selected points/particles
     uint32 size()const override
     {
         return selectedPoints_.size();
     }
 
+    /// @brief  Checks if there are no selected points
+    /// @return True if no points are selected, false otherwise
     bool empty()const override
     {
         return selectedPoints_.empty();
     }
 
+    /// @brief  Returns the volumes of the selected points (this is normally not used)
     span<const real> volumes()const override
     {
         return span<const real>(volume_.data(), volume_.size());
     }
 
-    span<const real> eqDiameters()const override
+    /// @brief  Returns the equivalent diameters of the regions (this is normally not used )
+        span<const real> eqDiameters()const override
     {
         return span<const real>(diameter_.data(), diameter_.size());
     }
 
+    /// @brief  Returns the center positions of the selected points
+    /// @return Span containing the center positions of all selected points
     span<const realx3> centers()const override
     {
         return span<const realx3>(center_.data(), center_.size());
     }
 
+    /// @brief  Returns the indices of the selected points (const version)
+    /// @param elem Element index (not used in this implementation)
+    /// @return Span containing the indices of all selected points
     span<const uint32> indices(uint32 elem)const override
     {
         return span<const uint32>(selectedPoints_.data(), selectedPoints_.size());
     }
 
+    /// @brief  Returns the indices of the selected points (non-const version)
+    /// @param elem Element index (not used in this implementation)
+    /// @return Span containing the indices of all selected points
     span<uint32> indices(uint32 elem) override
     {
         return span<uint32>(selectedPoints_.data(), selectedPoints_.size());
     }
 
-    /// @brief  update the selected points based on the ids
+    /// @brief  Updates the selected points based on the particle IDs
-    /// @return true if the operation is successful
+    /// @return True if the operation is successful, false otherwise
     bool update() override;
 
+    /// @brief  Checks if the data should be written to the same time file
+    /// @return True if data should be written to the same time file, false otherwise
     bool writeToSameTimeFile()const override
     {
         return true;
     }
 
+    /// @brief  Writes the data to the output stream
+    /// @param os Output stream
+    /// @return True if the operation is successful, false otherwise
     bool write(iOstream& os)const override;
 
 }; // class centerPointsRegionPoints

src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.cpp

@@ -30,7 +30,8 @@ pFlow::lineRegionPoints::lineRegionPoints
     if(raddi.size() != nPoints)
     {
         fatalErrorInFunction
-            << "The number elements in of radii list should be equal to the number of points"<<endl;
+            << "The number elements in of radii list should be equal to the "
+            << "number of points"<<endl;
             fatalExit;
     }
     
@@ -50,7 +51,7 @@ pFlow::lineRegionPoints::lineRegionPoints
 
 pFlow::span<const pFlow::uint32> pFlow::lineRegionPoints::indices(uint32 elem) const
 {
-    if(elem>=size())
+    if(elem >= size())
     {
         fatalErrorInFunction
             << "The element index is out of range. elem: " << elem
@@ -58,13 +59,15 @@ pFlow::span<const pFlow::uint32> pFlow::lineRegionPoints::indices(uint32 elem) c
         fatalExit;
     }
     
-    return span<const uint32>(selectedPoints_[elem].data(), selectedPoints_[elem].size());
+    return span<const uint32>(
+        selectedPoints_[elem].data(), 
+        selectedPoints_[elem].size());
 }
 
 bool pFlow::lineRegionPoints::update()
 {
     const auto points = database().updatePoints();
-    for(auto& elem:selectedPoints_)
+    for(auto& elem : selectedPoints_)
     {
         elem.clear();
     }
@@ -84,17 +87,18 @@ bool pFlow::lineRegionPoints::update()
 
 bool pFlow::lineRegionPoints::write(iOstream &os) const
 {
-    os <<"# Spheres along a straight line \n";
+    os << "# Spheres along a straight line \n";
-    os <<"#  No."<<tab <<"centerPoint" << tab <<"diameter"<<endl;
+    os << "#  No." << tab << "centerPoint" << tab << "diameter" << endl;
-    for(uint32 i=0; i< sphereRegions_.size(); ++i)
+    for(uint32 i=0; i < sphereRegions_.size(); ++i)
     {
-        os <<"# "<<i<<tab<<sphereRegions_[i].center() << tab <<diameters_[i] << '\n';
+        os << "# " << i << tab << sphereRegions_[i].center() 
+           << tab << diameters_[i] << '\n';
     }
-    os<<"time/No. ";
+    os << "time/No. ";
-    for(uint32 i=0; i< sphereRegions_.size(); ++i)
+    for(uint32 i=0; i < sphereRegions_.size(); ++i)
     {
-        os <<i<<tab;
+        os << i << tab;
     }
-    os <<endl;
+    os << endl;
     return true;
 }

src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.hpp

@@ -18,14 +18,42 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+/**
+ * @class lineRegionPoints 
+ * @brief Spherical regions along a line for selecting points/particles
+ * 
+ * The lineRegionPoints class is responsible for selecting points/particles along a 
+ * specified line and creating sphere regions around those points. It partitions
+ * the line into multiple sphere regions (equally spaced) and keeps track of 
+ * which points/particles in the simulation fall into each region.
+ * 
+ * This class is used for post-processing data by analyzing distributions of
+ * particles along a linear path through the simulation domain. It maintains:
+ *   - A line defining the sampling path
+ *   - Spherical regions along this line that include particles
+ *   - Center points for each region
+ *   - Volumes and diameters of regions
+ *   - Indices of points/particles contained in each region
+ * 
+ * The regions can be updated as the simulation progresses, and the data
+ * can be written to output for analysis.
+ *
+ * @see regionPoints
+ * @see line
+ * @see sphere
+ * @see fieldsDataBase
+ */
+
 #ifndef __lineRegionPoints_hpp__
 #define __lineRegionPoints_hpp__
 
+
 #include "regionPoints.hpp"
 #include "line.hpp"
 #include "sphere.hpp"
 #include "Vectors.hpp"
 
+
 namespace pFlow
 {
 
@@ -35,68 +63,80 @@ class lineRegionPoints
 {
 private:
 
-    /// line region for selecting points
+    /// Line path defining the axis of the spherical regions
     line                line_;
 
-    /// all sphere regions 
+    /// Collection of sphere regions along the line
     Vector<sphere>        sphereRegions_;
 
-    /// center poitns of regions/elements
+    /// Center points of all spherical regions
     realx3Vector          centerPoints_;
 
-    /// volumes of all elements/regions
+    /// Volumes of all spherical regions
     realVector            volumes_;
 
+    /// Diameters of all spherical regions
     realVector            diameters_;
 
-    /// the point indices that are selected by this region
+    /// Point/particles indices selected by each region
     Vector<uint32Vector>  selectedPoints_;
 
 public:
 
+    /// Type information for runtime type identification
     TypeInfo(line::TYPENAME());
 
+    /// Construct from dictionary that contains lineInfo and fields database
     lineRegionPoints(
         const dictionary& dict,
         fieldsDataBase& fieldsDataBase);
 
+    /// Default destructor
     ~lineRegionPoints() override = default;
 
+    /// Return number of regions
     uint32 size()const override
     {
         return sphereRegions_.size();
     }
 
+    /// Check if regions list is empty
     bool empty()const override
     {
         return sphereRegions_.empty();
     }
 
+    /// Return volumes of all regions
     span<const real> volumes()const override
     {
         return span<const real>(volumes_.data(), volumes_.size());
     }
 
+    /// Return equivalent diameters of all regions
     span<const real> eqDiameters()const override
     {
         return span<const real>(diameters_.data(), diameters_.size());
     }
 
+    /// Return center points of all regions
     span<const realx3> centers()const override
     {
         return span<const realx3>(centerPoints_.data(), centerPoints_.size());
     }
 
+    /// Return indices of points in the specified element/region
     span<const uint32> indices(uint32 elem)const override;
     
-
+    /// Update regions based on current particle positions
     bool update() override;
 
+    /// Whether to write all data to the same time file
     bool writeToSameTimeFile()const override
     {
         return true;
     }
 
+    /// Write data to output stream
     bool write(iOstream& os) const override;
 };
     

src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp

@@ -0,0 +1,98 @@
+#include "multipleSpheresRegionPoints.hpp"
+#include "fieldsDataBase.hpp"
+
+pFlow::multipleSpheresRegionPoints::multipleSpheresRegionPoints
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDataBase
+)
+:
+    regionPoints(dict, fieldsDataBase)
+{
+    const auto& multiSphereInfo = dict.subDict("multipleSphereInfo");
+    
+    // Read centers and radii lists
+    auto centers = multiSphereInfo.getVal<List<realx3>>("centers");
+    auto radii = multiSphereInfo.getVal<List<real>>("radii");
+    
+    // Check if lists have the same length
+    if(centers.size() != radii.size())
+    {
+        fatalErrorInFunction
+            << "The number of centers (" << centers.size() 
+            << ") does not match the number of radii (" << radii.size() << ")"
+            << endl;
+        fatalExit;
+    }
+    
+    uint32 nSpheres = centers.size();
+    
+    // Initialize data structures
+    sphereRegions_.resize(nSpheres, sphere(realx3(0.0, 0.0, 0.0), 1.0));
+    centerPoints_.resize(nSpheres);
+    diameters_.resize(nSpheres);
+    volumes_.resize(nSpheres);
+    selectedPoints_.resize(nSpheres);
+
+    // Setup each sphere
+    for (uint32 i = 0; i < nSpheres; ++i)
+    {
+        real diameter = 2.0 * radii[i]; // Convert radius to diameter
+        
+        sphereRegions_[i] = pFlow::sphere(centers[i], radii[i]);
+        centerPoints_[i] = centers[i];
+        diameters_[i] = diameter;
+        volumes_[i] = sphereRegions_[i].volume();
+    }
+}
+
+pFlow::span<const pFlow::uint32> pFlow::multipleSpheresRegionPoints::indices(uint32 elem) const
+{
+    if (elem >= size())
+    {
+        fatalErrorInFunction
+            << "The element index is out of range. elem: " << elem
+            << " size: " << size() << endl;
+        fatalExit;
+    }
+
+    return span<const uint32>(selectedPoints_[elem].data(), selectedPoints_[elem].size());
+}
+
+bool pFlow::multipleSpheresRegionPoints::update()
+{
+    const auto points = database().updatePoints();
+    for (auto& elem : selectedPoints_)
+    {
+        elem.clear();
+    }
+
+    for (uint32 i = 0; i < points.size(); ++i)
+    {
+        for (uint32 j = 0; j < sphereRegions_.size(); ++j)
+        {
+            if (sphereRegions_[j].isInside(points[i]))
+            {
+                selectedPoints_[j].push_back(i);
+            }
+        }
+    }
+    return true;
+}
+
+bool pFlow::multipleSpheresRegionPoints::write(iOstream &os) const
+{
+    os << "# Multiple spheres region points\n";
+    os << "#  No." << tab << "centerPoint" << tab << "diameter" << endl;
+    for (uint32 i = 0; i < sphereRegions_.size(); ++i)
+    {
+        os << "# " << i << tab << sphereRegions_[i].center() << tab << diameters_[i] << '\n';
+    }
+    os << "time/No. ";
+    for (uint32 i = 0; i < sphereRegions_.size(); ++i)
+    {
+        os << i << "  ";
+    }
+    os << endl;
+    return true;
+}

src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.hpp

@@ -0,0 +1,161 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+
+/**
+ * @class       multipleSpheresRegionPoints
+ * @brief       A class to select and track particles contained within multiple 
+ *              spherical regions
+ * @details     This class defines multiple spherical regions in the simulation 
+ *             domain and identifies which particles are contained within each 
+ *             sphere at each time step. It inherits from the regionPoints base 
+ *             class and specializes it for handling multiple spherical regions 
+ *             simultaneously.
+ *             
+ *             The class reads a list of center points and radii from a dictionary,
+ *             creates sphere objects for each, and provides methods to:
+ *             - Track which particles are inside each spherical region
+ *             - Return volumetric information about the regions
+ *             - Access center points and size information for each sphere
+ *
+ *             
+ * @note        Used in post-processing workflows to analyze particle behavior
+ *             in specific regions of interest within the simulation domain.
+ * 
+ * @see        regionPoints Base class for all region-based point selection
+ * @see        sphere Geometric primitive used to define spherical regions
+ * @see        postprocessPhasicFlow Utility for post-processing simulation data
+ * @see        fieldsDataBase Database containing simulation field data
+ */
+
+#ifndef __multipleSpheresRegionPoints_hpp__
+#define __multipleSpheresRegionPoints_hpp__
+
+#include "regionPoints.hpp"
+#include "sphere.hpp"
+#include "Vectors.hpp"
+
+namespace pFlow
+{
+
+class multipleSpheresRegionPoints
+:
+    public regionPoints
+{
+private:
+
+    /// Vector containing all spherical regions used for particle selection
+    Vector<sphere>        sphereRegions_;
+
+    /// Center coordinates of all spherical regions
+    realx3Vector          centerPoints_;
+
+    /// Diameters of all spherical regions
+    realVector            diameters_;
+
+    /// Volumes of all spherical regions
+    realVector            volumes_;
+
+    /// Vectors of point indices for particles contained in each spherical region
+    /// Each element corresponds to a particular sphere region
+    Vector<uint32Vector>  selectedPoints_;
+
+public:
+
+    /// Type identification for run-time type information
+    TypeInfo("multipleSpheres");
+
+    /// Constructor
+    /// @param dict Dictionary containing multipleSpheresInfo for the regions
+    /// @param fieldsDataBase Reference to the database containing field data
+    multipleSpheresRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    /// Virtual destructor for proper inheritance cleanup
+    ~multipleSpheresRegionPoints() override = default;
+
+    /// Returns the number of spherical regions
+    /// @return Number of spherical regions
+    uint32 size()const override
+    {
+        return sphereRegions_.size();
+    }
+
+    /// Checks if there are any spherical regions defined
+    /// @return True if no regions exist, false otherwise
+    bool empty()const override
+    {
+        return sphereRegions_.empty();
+    }
+
+    /// Returns the volumes of all spherical regions
+    /// @return Span containing the volumes of all regions
+    span<const real> volumes()const override
+    {
+        return span<const real>(volumes_.data(), volumes_.size());
+    }
+
+    /// Returns the diameters of all spherical regions
+    /// @return Span containing the diameters of all regions
+    span<const real> diameters()const
+    {
+        return span<const real>(diameters_.data(), diameters_.size());
+    }
+
+    /// Returns the equivalent diameters of all spherical regions
+    /// @return Span containing the equivalent diameters (same as diameters)
+    span<const real> eqDiameters()const
+    {
+        return diameters();
+    }
+
+    /// Returns the center coordinates of all spherical regions
+    /// @return Span containing the center points of all regions
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(centerPoints_.data(), centerPoints_.size());
+    }
+
+    /// Returns the indices of particles contained in a specific spherical region
+    /// @param elem Index of the spherical region to query
+    /// @return Span containing indices of particles within the specified region
+    span<const uint32> indices(uint32 elem)const override;
+
+    /// Updates the selection of particles within each spherical region
+    /// @return True if update was successful, false otherwise
+    bool update() override;
+
+    /// Determines if data should be written to the same time file
+    /// @return True to indicate regions should be written to the same time file
+    bool writeToSameTimeFile()const override
+    {
+        return true;
+    }
+
+    /// Writes region data to the output stream
+    /// @param os Output stream to write data to
+    /// @return True if write operation was successful, false otherwise
+    bool write(iOstream& os) const override;
+};
+    
+}
+
+
+#endif // __multipleSpheresRegionPoints_hpp__

src/PostprocessData/region/regionPoints/regionPoints/regionPoints.hpp

@@ -32,27 +32,43 @@ namespace pFlow
 class fieldsDataBase;
 class Time;
 
+/**
+ * @class regionPoints
+ * @brief Abstract base class for managing and processing volumetric regions 
+ * in the simulation. Particles are selected based on their positions within
+ * these defined regions.
+ * 
+ * This class provides an interface for accessing and manipulating data 
+ * related to regions of points (particles), including their volumes, equivalent 
+ * diameters, center points, and particle indices that they contain. It interacts with the 
+ * fieldsDataBase and Time classes to retrieve simulation data. 
+ */
 class regionPoints
 {
-    using PointsTypeHost = typename pointStructure::PointsTypeHost;
+        using PointsTypeHost = typename pointStructure::PointsTypeHost;
 
+    /// Reference to the fields database containing simulation data
     fieldsDataBase& fieldsDataBase_;
 
 public:
 
     TypeInfo("regionPoints");
 
-
+    /// Constructor with dictionary and fields database reference
     regionPoints(
         const dictionary& dict,
         fieldsDataBase& fieldsDataBase);
     
+    /// Default virtual destructor
     virtual ~regionPoints() = default;
 
+    /// Returns reference to the time object from the database
     const Time& time()const;
 
+    /// Returns const reference to the fields database
     const fieldsDataBase& database()const;
 
+    /// Returns non-const reference to the fields database
     fieldsDataBase& database();
     
     /// @brief size of elements 
@@ -61,11 +77,7 @@ public:
     
     /// @brief  check if the region is empty
     virtual 
-    bool empty()const = 0;
+    bool empty()const = 0;  
-
-    /*/// @brief return the type of the region
-    virtual const word& type()const = 0;*/
-    
     
     /// @brief  volume of elements
     /// @return sapn for accessing the volume of elements 
@@ -75,30 +87,29 @@ public:
     virtual
     span<const real> eqDiameters()const = 0; 
 
-    /// center points of elements 
+    /// center points of elements
     virtual
     span<const realx3> centers()const = 0;
 
-    /// indices of particles inside the element @var elem
+    /// Returns const span of particle indices inside the specified element region
     virtual 
     span<const uint32> indices(uint32 elem)const = 0;
 
+    /// Returns non-const span of particle indices inside the specified element region
     virtual
     span<uint32> indices(uint32 elem) = 0;
 
+    /// Updates the points (particles) inside regions based on current particle positions
     virtual 
     bool update() = 0;
     
+    /// Returns true if the region should be written to the same time file
     virtual
     bool writeToSameTimeFile()const = 0;
 
+    /// Writes region data to the output stream
     virtual 
     bool write(iOstream& os)const=0;
-
-    /*static
-    uniquePtr<regionPoints> create(
-        const dictionary& dict,
-        fieldsDataBase& fieldsDataBase);*/
     
 };
 

src/PostprocessData/region/regionPoints/sphereRegionPoints/sphereRegionPoints.hpp

@@ -18,6 +18,19 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+/**
+ * @file sphereRegionPoints.hpp
+ * @brief A class representing a spherical region for point selection
+ * 
+ * This class provides functionality to select points within a spherical region
+ * and to compute related properties such as volume and equivalent diameter.
+ * It inherits from regionPoints and implements all required virtual methods.
+ * 
+ * @see regionPoints
+ * @see sphere
+ * @see fieldsDataBase
+ */
+
 #ifndef __sphereRegionPoints_hpp__
 #define __sphereRegionPoints_hpp__
 
@@ -27,75 +40,127 @@ Licence:
 
 namespace pFlow
 {
+
 class sphereRegionPoints
 :
     public regionPoints
 {
 private:
 
-    /// spehre region for selecting points
+    /// Sphere object defining the region for point selection
     sphere  sphereRegion_;
 
-    /// the volume of region
+    /// Volume of the spherical region
     real    volume_;
 
+    /// Diameter of the spherical region
     real    diameter_;
 
-    /// the point indices that are selected by this region
+    /// Indices of points that are selected by this region
     uint32Vector    selectedPoints_;
 
 public:
 
     TypeInfo(sphere::TYPENAME());
 
+    /**
+     * @brief Construct a spherical region for point selection
+     * 
+     * @param dict Dictionary containing sphereInfo dictionary 
+     * @param fieldsDataBase Database containing fields data
+     */
     sphereRegionPoints(
         const dictionary& dict,
         fieldsDataBase& fieldsDataBase);
 
+    /// Destructor
     ~sphereRegionPoints() override = default;
     
+    /**
+     * @brief Get the number of regions (always 1 for sphere)
+     * @return Always returns 1
+     */
     uint32 size()const override
     {
         return 1;
     }
 
+    /**
+     * @brief Check if the region is empty
+     * @return Always returns false
+     */
     bool empty()const override
     {
         return false;
     }
 
+    /**
+     * @brief Get the volume of the spherical region
+     * @return A span containing the volume of the region
+     */
     span<const real> volumes()const override
     {
         return span<const real>(&volume_, 1);
     }
 
+    /**
+     * @brief Get the equivalent diameter of the spherical region
+     * @return A span containing the diameter of the region
+     */
     span<const real> eqDiameters()const override
     {
         return span<const real>(&diameter_, 1);
     }
 
+    /**
+     * @brief Get the center of the spherical region
+     * @return A span containing the center point of the region
+     */
     span<const realx3> centers()const override
     {
         return span<const realx3>(&sphereRegion_.center(), 1);
     }
-        
+    
+    /**
+     * @brief Get the indices of points within the region (const version)
+     * @param elem Element index (ignored as there's only one sphere)
+     * @return A span containing indices of points within the region
+     */    
     span<const uint32> indices(uint32 elem)const override
     {
         return span<const uint32>(selectedPoints_.data(), selectedPoints_.size());
     }
 
+    /**
+     * @brief Get the indices of points within the region (non-const version)
+     * @param elem Element index (ignored as there's only one sphere)
+     * @return A span containing indices of points within the region
+     */
     span<uint32> indices(uint32 elem) override
     {
         return span<uint32>(selectedPoints_.data(), selectedPoints_.size());
     }
 
+    /**
+     * @brief Update the points selected by this region
+     * @return True if update was successful
+     */
     bool update()override;
 
+    /**
+     * @brief Determine if data should be written to the same time file
+     * @return Always returns true
+     */
     bool writeToSameTimeFile()const override
     {
         return true;
     }
 
+    /**
+     * @brief Write region data to output stream
+     * @param os Output stream to write to
+     * @return True if write was successful
+     */
     bool write(iOstream& os)const override;
     
 };

src/PostprocessData/sampleDictionary/postprocessDataDict

@@ -2,107 +2,128 @@
 |  phasicFlow File                                                            | 
 |  copyright: www.cemf.ir                                                     | 
 \* ------------------------------------------------------------------------- */  
-objectName 	processDataDict;
+objectName   postprocessDataDict;
-objectType 	dictionary;;
+objectType   dictionary;;
-fileFormat  ASCII;
+fileFormat   ASCII;
 /*---------------------------------------------------------------------------*/
 
-runTimeActive 		yes;
+runTimeActive yes;
 
 defaultTimeControl
 {
-	timeControl;
+	timeControl         timeStep;
-	startTime;
+	startTime           0;
-	endTime;
+	endTime             1000;
-	actionInterval  0.05; 
+	executionInterval   150; 
 }
 
 components
 (
-    velocityProb
+      velocityProb
-    {
+      {
-        method 	particleProbe;
+          processMethod    particleProbe;
-        region  idSelecttion;
+          processRegion    centerPoints;
-        field   velocity;
+          selector        id;
-        ids (0 10 100);
+          field   component(position,y);
-        timeControl 	timeStep;
+          ids (0 10 100);
-        startTime 		0;
+      }
-        endTime 		infinity;
+    
-        probInterval 	1;
+      onSingleSphere
-    }
+      {
+          // method of performing the sum (arithmetic, uniformDistribution, GaussianDistribution)
+          processMethod    arithmetic; 
+          processRegion    sphere; // type of region on which processing is performed 
+          
+          sphereInfo
+          {
+              radius 0.01;
+              center (-0.08 -0.08 0.015);
+          }
+          
+          timeControl default; // settings, timeStep, simulationTime 
+          
+          /// all the post process operations to be done 
+          operations
+          (
+              // computes the arithmetic mean of particle velocity 
+              averageVel
+              {
+                  function	  average;
+                  field           velocity;
+                  dividedByVolume 	  no;  //default 
+                  threshold 	  3;   //default is 1;
+                  includeMask     all;
+              }
+              
+              // computes the fraction of par1 in the region 
+              par1Fraction
+              {
+                  function        average;
+                  field           one;
+                  phi             one; // default
+                  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;  
+                  }
+              }
 
-    comp2
+              numberDensity
-    {
+              {
-        method    uniformDistribution;
+                  function        sum;
-        region spehre;
+                  field           one;
-        
+                  phi             one; // default
-        sphereInfo
+                  dividedByVolume yes;
-        {
+                  
-            radius 0.01;
+              }
-            center ();
+          );
-        }
-        timeControl default; //default;        
-        operations
-        (
-            numParticle
-            {
-                function		sum;
-                field           compoenent(velocity,x);
-                phi             square(mass);
-                divideByVol 	no;  //default 
-                threshold 		1;   //default;
-                defaultVal      NaN;
-                //includeMask 	all; //default;
-                includeMask     lessThan;
-                lessThanInfo
-                {
-                    field      diameter;
-                    value      0.003;
-                }
-            }
-        );
     }
-
+    
-    comp3
+    alongALine
     {
-        
+          processMethod    arithmetic; 
-        region  line;
+          processRegion    line;
-        lineInfo
+          
-        {
+          // the time interval for executing the post-processing 
-            p1 ();
+          // other options: timeStep, default, and settings
-            p2 ();
+          timeControl simulationTime; 
+          startTime   1.0;
+          endTime     3.0;
+          executionInterval 0.1;
+           
+          
+          // 10 spheres with radius 0.01 along the straight line defined by p1 and  p2 
+          lineInfo
+          {
+            p1 (0 0 0);
+            p2 (0 0.15 0.15);
             numPoints 10;
-            radius 0.01;
+            radius    0.01;
-        }
+          }
-        timeControl settingsDict; //default; 
+          
-        type    numberBased;
+          operations
-        operations();
+          (
+              // computes the arithmetic mean of particle velocity 
+              numberDensity
+              {
+                  function	   sum;
+                  field            one;
+                  dividedByVolume  yes;  //default is no 
+              }
+              
+              volumeDensity
+              {
+                function         sum;
+                field            cube(diameter); // d^3, although it differs by pi/6
+                dividedByVolume  yes;  //default is no 
+              }
+         );
     }
-
+    
-    comp4
+);
-    {
-        type    GaussianDistribution;
-        region 	hexMesh; // unstructuredMehs; 
-        hexMeshInfo
-        {
-            min (-0.3 -1.4 -0.01);
-            max ( 0.3  2  0.48  );
-            nx 30; 		// number of divisions in x direction
-            ny 160; 	// number of divisions in y direction
-            nz 24; 		// number of divisions in z direction
-        }
-        timeControl settingsDict; // read from settingsDict
-        operations
-        (
-            avVelocity
-            {
-                type 		average;
-                field  		realx3 velocity; // default to real 1.0 
-                divideByVol no; // default
-                threshold   1;  //default;
-                includeMask all; //default;
-                
-            }
-        );
-    }
-);