first commit for post-processing

- the whole structure is ready.
- next step whould be execute methods and then write methods
- post-processing after simulation is not started yet.

src/PostProcessData/CMakeLists.txt

@@ -0,0 +1,35 @@
+
+set(SourceFiles 
+fieldsDataBase/fieldsDataBase.cpp
+
+postprocessComponent/postprocessComponent.cpp
+postprocessComponent/PostprocessComponents.cpp
+
+operation/postprocessOperation.cpp
+operation/PostprocessOperationSum.cpp
+
+
+postprocessData.cpp
+
+postprocessComponent/postprocessComponent.cpp
+postprocessComponent/PostprocessComponents.cpp
+postprocessComponent/particleProbePostprocessComponent.cpp 
+
+operation/includeMask/includeMask.cpp
+operation/includeMask/IncludeMasks.cpp
+operation/postprocessOperation.cpp
+operation/PostprocessOperationSum.cpp
+
+region/regionPoints/sphereRegionPoints.cpp
+region/regionPoints/regionPoints.cpp
+region/regionPoints/lineRegionPoints.cpp
+region/regionPoints/centerPointsRegionPoints.cpp
+
+)
+
+set(link_libs Kokkos::kokkos phasicFlow Particles)
+
+pFlow_add_library_install(PostProcessData SourceFiles link_libs)
+
+add_subdirectory(testPostprocess)
+

src/PostProcessData/fieldsDataBase/fieldFunctions.hpp

@@ -0,0 +1,165 @@
+/*------------------------------- 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 __fieldFunctions_hpp__
+#define __fieldFunctions_hpp__
+
+#include "types.hpp"
+#include "span.hpp"
+
+namespace pFlow
+{
+
+template<typename T>
+inline
+void funcCast(span<T> src, span<real> dst)
+{
+    for(uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = static_cast<real>(src[i]);
+    }
+}
+
+template<typename T>
+inline
+void funcAbs(span<T> src, span<real> dst)
+{
+    for(uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = abs(src[i]);
+    }
+}
+
+template<typename T>
+inline
+void funcSquare(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = pow(static_cast<real>(src[i]),2);
+    }
+}
+
+template<typename T>
+inline
+void funcCube(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = pow(static_cast<real>(src[i]),3);
+    }
+}
+
+template<typename T>
+inline
+void funcSquareRoot(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = sqrt(static_cast<real>(src[i]));
+    }
+}
+
+template<VectorType T> 
+inline
+void funcMagnitude(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = src[i].length();
+    }
+}
+
+template<VectorType T>
+inline
+void funcMagnitudeSquare(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = dot(src[i],src[i]);
+    }
+}
+
+template<VectorType T>
+inline
+void funcMagnitudeCube(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = pow(src[i].length(),3);
+    }
+}
+
+template<VectorType T>
+inline 
+void funcMagnitudeSquareRoot(span<T> src, span<real> dst)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        dst[i] = sqrt(src[i].length());
+    }
+}
+
+inline
+void funcComponent(span<realx3> src, span<real> dst, char component)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        switch(component)
+        {
+            case 'x':
+                dst[i] = src[i].x();
+                break;
+            case 'y':
+                dst[i] = src[i].y();
+                break;
+            case 'z':
+                dst[i] = src[i].z();
+                break;
+        }
+    }
+}
+
+inline
+void funcComponent(span<realx4> src, span<real> dst, char component)
+{
+    for( uint32 i=0; i<src.size(); i++)
+    {
+        switch(component)
+        {
+            case 'x':
+                dst[i] = src[i].x();
+                break;
+            case 'y':
+                dst[i] = src[i].y();
+                break;
+            case 'z':
+                dst[i] = src[i].z();
+                break;
+            case 'w':
+                dst[i] = src[i].w();
+                break;
+        }
+    }
+}
+
+} // namespace pFlow
+
+#endif //__fieldFunctions_hpp__

src/PostProcessData/fieldsDataBase/fieldsDataBase.cpp

@@ -0,0 +1,664 @@
+/*------------------------------- 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 <regex>
+
+#include "vocabs.hpp"
+#include "Time.hpp"
+#include "fieldsDataBase.hpp"
+#include "fieldFunctions.hpp"
+#include "dynamicPointStructure.hpp"
+
+bool pFlow::fieldsDataBase::checkForUpdate(const word &name, bool forceUpdate)
+{
+    auto t = currentTime();
+    bool shouldUpdate = false;
+    
+    if(auto [iter, found]= captureTime_.findIf(name); found)
+    {
+        shouldUpdate = iter->second < t;
+        iter->second = t;
+    }
+    else
+    {
+        shouldUpdate = true;
+        captureTime_.insertIf(name, t);
+    }
+
+    return shouldUpdate;
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetRealField(const word &name)
+{
+    
+    bool shouldUpdate = checkForUpdate(name);
+        
+    if(shouldUpdate)
+    {
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            name,
+            FieldTypeHost<real>
+            (
+                name, 
+                "value",
+                pointFieldSize()
+            )
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(name);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+bool pFlow::fieldsDataBase::findFunction(
+    const word &compoundFieldName,
+    word &fieldName,
+    fieldsDataBase::Functions &func)
+{
+
+    std::regex pattern(R"((\w+)?\((\w+)(?:,([xyzw]))?\)|(\w+))");
+    std::smatch match;
+
+    if (std::regex_match(compoundFieldName, match, pattern))
+    {
+        if (match[1].matched) // Function is present
+        {
+            word functionName = match[1].str();
+            fieldName = match[2].str();
+
+            // Map the function name to the enum value
+            if(functionName=="component")
+            {
+                if (!match[3].matched) // Component is not present
+                {
+                    fatalErrorInFunction << 
+                        "Component (x, y, z, or w) is not specified in the function component: " << compoundFieldName <<
+                        " the format is component(u,x), where u is the vector field name and x is the compoenent." << endl;
+                    return false;
+                }
+
+                switch (match[3].str()[0])
+                {
+                    case 'x': func = fieldsDataBase::Functions::ComponentX; break;
+                    case 'y': func = fieldsDataBase::Functions::ComponentY; break;
+                    case 'z': func = fieldsDataBase::Functions::ComponentZ; break;
+                    case 'w': func = fieldsDataBase::Functions::ComponentW; break;
+                    default:
+                        fatalErrorInFunction << 
+                            "Invalid component specified: " << match[3].str() << endl;
+                        return false;
+                }
+                return true;       
+            }
+            else if (functionName == "abs")
+            {
+                func = fieldsDataBase::Functions::Abs;
+            }
+            else if (functionName == "square")
+            {
+                func = fieldsDataBase::Functions::Square;
+            }
+            else if (functionName == "cube")
+            {
+                func = fieldsDataBase::Functions::Cube;
+            }
+            else if (functionName == "sqrt")
+            {
+                func = fieldsDataBase::Functions::SqureRoot;
+            }
+            else if (functionName == "mag")
+            {
+                func = fieldsDataBase::Functions::Magnitude;
+            }
+            else if (functionName == "magSquare")
+            {
+                func = fieldsDataBase::Functions::MagnitudeSquare;
+            }
+            else if (functionName == "magCube")
+            {
+                func = fieldsDataBase::Functions::MagnitudeCube;
+            }
+            else if (functionName == "magSqrt")
+            {
+                func = fieldsDataBase::Functions::MagnitudeSquareRoot;
+            }
+            else
+            {
+                fatalErrorInFunction << 
+                    "Unknown function specified: " << functionName<< 
+                    " in: "<<compoundFieldName << endl;
+                return false;
+            }
+
+            if(match[3].matched)
+            {
+                fatalErrorInFunction << 
+                    "The function: " << functionName << 
+                    " does not accept component (x, y, z, or w) as the second argument (in: "
+                    << compoundFieldName<<" )." << endl;
+                return false;
+            }
+            else
+            {
+                return true;    
+            }
+            
+        }
+        else if (match[4].matched) // Only fieldName is present
+        {
+            fieldName = match[4].str();
+            func = fieldsDataBase::Functions::None; // No function
+            return true;
+        }
+        
+    }
+
+    return false; // No match
+}
+
+bool pFlow::fieldsDataBase::inputOutputType
+(
+    fieldsDataBase::Functions func, 
+    const word &inputType, 
+    word &outputType
+)
+{
+    switch (func)
+    {
+        case Functions::ComponentX:
+        case Functions::ComponentY:
+        case Functions::ComponentZ:
+            if(inputType == getTypeName<realx3>() ||
+               inputType == getTypeName<realx4>() )
+            {
+                outputType = getTypeName<real>();
+                return true;
+            }
+            else
+            {
+                fatalErrorInFunction<<"Wrong function: component(u,comp), for input field type: "<<inputType<<endl; 
+                return false;
+            }
+        case Functions::ComponentW:
+            if(inputType == getTypeName<realx4>() )
+            {
+                outputType = getTypeName<real>();
+                return true;
+            }
+            else
+            {
+                fatalErrorInFunction<<"Wrong function: component(u,w), for input field type: "<<inputType<<endl; 
+                return false;
+            }
+            break;
+        case Functions::Abs:
+        case Functions::Square:
+        case Functions::Cube:
+        case Functions::SqureRoot:
+            if(inputType == getTypeName<real>() || 
+               inputType == getTypeName<uint32>() ||
+               inputType == getTypeName<int32>() ||
+               inputType == getTypeName<int64>() )
+            {
+                outputType = getTypeName<real>();
+                return true;
+            }
+            else
+            {
+                fatalErrorInFunction<<"Wrong input field type for functions abs, squqre, cube, and sqrt."<<
+                    " field type is "<< inputType<<endl;
+                return false;
+            }
+            break;
+        case Functions::Magnitude:
+        case Functions::MagnitudeSquare:
+        case Functions::MagnitudeCube:
+        case Functions::MagnitudeSquareRoot:
+            if(inputType == getTypeName<realx3>() ||
+               inputType == getTypeName<realx4>() )
+            {
+                outputType = getTypeName<real>();
+                return true;
+            }
+            else
+            {
+                fatalErrorInFunction<<"Wroing input field type for functions mag, magSquare, magCube, magSqrt. "<<
+                    " Input field type is "<< inputType<<endl;
+                return false;
+            }
+            break;
+        case Functions::None:
+            if(inputType == getTypeName<realx3>() ||
+               inputType == getTypeName<realx4>() ||
+               inputType == getTypeName<real>())
+            {
+                outputType = inputType;
+                return true;
+            }
+            else if( inputType == getTypeName<uint32>() ||
+                     inputType == getTypeName<int32>() )
+            {
+                outputType = getTypeName<real>();
+                return true;
+            }
+            else
+            {
+                fatalErrorInFunction<< "Wroing input field type "<< inputType<<endl;
+                return false;
+            }
+            break;
+        default:
+            fatalErrorInFunction<<"Wroing function"<<endl;
+
+    }
+    return false;
+}
+
+pFlow::fieldsDataBase::fieldsDataBase(const Time &time, bool inSimulation)
+: 
+    time_(time), 
+    inSimulation_(inSimulation)
+{}
+
+pFlow::timeValue pFlow::fieldsDataBase::currentTime() const
+{
+    return time_.currentTime();
+}
+
+bool pFlow::fieldsDataBase::getPointFieldType
+(
+    const word& compoundName, 
+    word& fieldName,
+    word& originalType, 
+    word& typeAfterFunction,
+    Functions& func
+)
+{
+    if( !findFunction(compoundName, fieldName, func))
+    {
+        fatalErrorInFunction;
+        return false;
+    }
+
+    if(reservedFieldNames_.contains(fieldName))
+    {
+        originalType = reservedFieldNames_.find(fieldName)->second;
+    }
+    else
+    {
+        word fieldTypeName =  time_.lookupObjectTypeName(fieldName);
+        word space;
+        if(!pointFieldGetType(fieldTypeName, originalType, space))
+        {
+            fatalErrorInFunction<<"Cannot extract type name from "<< fieldTypeName<<endl;
+            return false;
+        }
+    }
+    
+    word outputType;
+    if(!inputOutputType(func, originalType, outputType))
+    {
+        fatalErrorInFunction<<"Cannnot determine the input and output types for: "<< compoundName<<endl;
+        return false;
+    }
+    typeAfterFunction = outputType;
+    return true;
+}
+
+bool pFlow::fieldsDataBase::getPointFieldType
+(
+    const word &compoundName, 
+    word &originalType, 
+    word &typeAfterFunction
+)
+{
+    Functions func;
+    word fieldName;
+    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+}
+
+bool pFlow::fieldsDataBase::getPointFieldType
+(
+    const word &compoundName, 
+    word &typeAfterFunction
+)
+{
+    Functions func;
+    word originalType;
+    word fieldName;
+    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+}
+
+pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updatePoints(bool forceUpdate)
+{
+    
+    bool shouldUpdate = checkForUpdate("position", forceUpdate);
+        
+    if(shouldUpdate)
+    {
+        const auto& pstruct = pStruct(); 
+        allFields_.emplaceBackOrReplace<PointsTypeHost>(
+            "position", 
+            pstruct.activePointsHost());
+    }
+   
+    auto& points = allFields_.getObject<PointsTypeHost>("position");
+    
+    return span<realx3>(
+        points.data(), 
+        points.size());
+
+}
+
+pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updateFieldRealx3
+(
+    const word &compoundName, 
+    bool forceUpdate
+)
+{
+    word originalType, typeAfterFunction, fieldName;
+    Functions func;
+    
+    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<realx3>(nullptr, 0);
+    }
+
+    if( originalType == getTypeName<realx3>() && func == Functions::None )
+    {
+        return updateField<realx3>(fieldName, forceUpdate);
+    }
+    else
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<realx3>(nullptr, 0);
+    }
+        
+}
+
+pFlow::span<pFlow::realx4> pFlow::fieldsDataBase::updateFieldRealx4
+(
+    const word &compoundName, 
+    bool forceUpdate
+)
+{
+    word originalType, typeAfterFunction, fieldName;
+    Functions func;
+    
+    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<realx4>(nullptr, 0);
+    }
+
+    if( originalType == getTypeName<realx4>() && func == Functions::None )
+    {
+        return updateField<realx4>(fieldName, forceUpdate);
+    }
+    else
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<realx4>(nullptr, 0);
+    }
+    
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::updateFieldReal
+(
+    const word &compoundName, 
+    bool forceUpdate
+)
+{
+    word originalType, typeAfterFunction, fieldName;
+    Functions func;
+    
+    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<real>(nullptr, 0);
+    }
+
+    // if the output type is not real, then it is not supported yet
+    if(typeAfterFunction != getTypeName<real>())
+    {
+        fatalErrorInFunction<< "The output type of field "<< compoundName<<
+            " is not real, it is: "<< typeAfterFunction<<endl;
+        fatalExit;
+        return span<real>(nullptr, 0);
+    }
+
+    // if the orginal type is real and no function, then it is a normal field
+    if( originalType == getTypeName<real>() && func == Functions::None )
+    {
+        return updateField<real>(fieldName, forceUpdate);
+    }
+    
+    // if the original type is uint32, and no funciton, cast to real 
+    if( originalType == getTypeName<uint32>() && func == Functions::None )
+    {
+        auto sp = updateField<uint32>(fieldName, forceUpdate);
+        auto spReal = createOrGetRealField(fieldName+".real");
+        funcCast(sp, spReal);
+        return spReal;
+    }
+    else
+    {
+        fatalErrorInFunction<<"No function can be applied to field of type uint32. "<<
+            " The field is: "<< compoundName<<endl;
+        return span<real>(nullptr, 0);
+        fatalExit;
+    }
+
+    if( originalType == getTypeName<real>() )
+    {
+        switch(func)
+        {
+            case Functions::Abs:
+            {
+                auto sp = updateField<real>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcAbs(sp, spReal);
+                return spReal;
+            }
+            case Functions::Square:
+            {
+                auto sp = updateField<real>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcSquare(sp, spReal);
+                return spReal;
+            }
+            case Functions::Cube:
+            {
+                auto sp = updateField<real>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcCube(sp, spReal);
+                return spReal;
+            }
+            case Functions::SqureRoot:
+            {
+                auto sp = updateField<real>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName+".sqrt");
+                funcSquareRoot(sp, spReal);
+                return spReal;
+            }
+            default:
+            {
+                fatalErrorInFunction<< "Wrong function for field type real in :"<<
+                    compoundName<<endl;
+                fatalExit;
+                return span<real>(nullptr, 0);
+            }
+        }
+    }
+
+    if( originalType == getTypeName<realx3>())
+    {
+        switch(func)
+        {
+            case Functions::Magnitude:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcMagnitude(sp, spReal);
+                return spReal;
+            }
+            case Functions::MagnitudeSquare:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcMagnitudeSquare(sp, spReal);
+                return spReal;
+            }
+            case Functions::MagnitudeCube:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcMagnitudeCube(sp, spReal);
+                return spReal;
+            }
+            case Functions::MagnitudeSquareRoot:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcMagnitudeSquareRoot(sp, spReal);
+                return spReal;
+            }
+            case Functions::ComponentX:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcComponent(sp, spReal, 'x');
+                return spReal;
+            }
+            case Functions::ComponentY:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcComponent(sp, spReal, 'y');
+                return spReal;
+            }
+            case Functions::ComponentZ:
+            {
+                auto sp = updateField<realx3>(fieldName, forceUpdate);
+                auto spReal = createOrGetRealField(compoundName);
+                funcComponent(sp, spReal, 'z');
+                return spReal;
+            }
+            default:
+            {
+                fatalErrorInFunction<< "Wrong function for field type realx3 in :"<<
+                    compoundName<<endl;
+                fatalExit;
+                return span<real>(nullptr, 0);
+            }
+        }
+    }
+
+    fatalErrorInFunction<<"NOT SUPPORTED "<<compoundName<<endl;
+    fatalExit;
+    // This line should never be reached
+    return span<real>(nullptr, 0);
+}
+
+pFlow::span<pFlow::uint32> pFlow::fieldsDataBase::updateFieldUint32
+(
+    const word& name,
+    bool forceUpdate
+)
+{
+    return updateField<uint32>(name, forceUpdate);
+}
+
+pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
+(
+    const word &compoundName, 
+    bool forceUpdate
+)
+{
+    word originalType, typeAfterFunction;
+    
+    if( !getPointFieldType(compoundName, originalType, typeAfterFunction))
+    {
+        fatalErrorInFunction<< "Error in getting the type name of field: "<<
+            compoundName<<", with type name: "<< originalType <<endl;
+        fatalExit;
+        return span<real>(nullptr,0);
+    }
+
+    if(  typeAfterFunction== getTypeName<realx3>()  )
+    {
+        return updateField<realx3>(compoundName, forceUpdate);
+    }
+    else if( typeAfterFunction == getTypeName<realx4>() )
+    {
+        return updateField<realx4>(compoundName, forceUpdate);
+    }
+    else if( typeAfterFunction == getTypeName<real>() )
+    {
+        return updateField<real>(compoundName, forceUpdate);
+    }
+    else
+    {
+        fatalErrorInFunction<< "Invalid feild "<< compoundName<<endl;
+        fatalExit;
+        return span<real>(nullptr, 0);
+    }
+}
+
+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)
+{
+	std::regex match("pointField\\<([A-Za-z1-9_]*)\\,([A-Za-z1-9_]*)\\>");
+    std::smatch search;
+	if(!std::regex_match(TYPENAME, search, match)) return false;
+    if(search.size()!= 3) return false;
+    fieldType  = search[1];
+    fieldSpace = search[2];
+	return true;
+}
+

src/PostProcessData/fieldsDataBase/fieldsDataBase.hpp

@@ -0,0 +1,389 @@
+/*------------------------------- 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 __fieldsDataBase_hpp__
+#define __fieldsDataBase_hpp__
+
+#include <variant>
+#include <concepts>
+
+
+#include "pointStructure.hpp"
+#include "pointFields.hpp"
+#include "anyList.hpp"
+#include "Map.hpp"
+#include "span.hpp"
+
+namespace pFlow
+{
+
+class Time;
+
+bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+
+template<typename T>
+concept ValidFieldType = 
+    std::same_as<T, real> || 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4> || 
+    std::same_as<T, uint32>;
+
+template<typename T> 
+concept VectorType = 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4>;
+
+template<typename T>
+concept ScalarType = 
+    std::same_as<T, real>;
+
+
+template<typename T>
+concept ValidRegionFieldType = 
+    std::same_as<T, real> || 
+    std::same_as<T, realx3> || 
+    std::same_as<T, realx4> ;
+
+
+using allPointFieldTypes = std::variant<span<real>, span<realx3>, span<realx4>>; 
+
+
+class fieldsDataBase
+
+{
+private:
+
+    // - Typedefs
+
+        /// Point position data type
+        using PointsTypeHost    = typename pointStructure::PointsTypeHost;
+        
+        /// Point field data type 
+        template<typename T>
+        using FieldTypeHost     = typename internalField<T>::FieldTypeHost;
+
+        /// @brief Enumeration of functions that can be applied to point fields.
+        enum class Functions
+        {
+            None,       /// no function
+            ComponentX, /// component(u,x)
+            ComponentY, /// component(u,y)
+            ComponentZ, /// component(u,z)
+            ComponentW, /// component(u,w)
+            Abs,        /// abs(s)
+            Square,     /// square(s)
+            Cube,       /// cube(s)
+            SqureRoot,  /// sqrt(s)
+            Magnitude,          /// mag(u)
+            MagnitudeSquare,    /// magSquare(u)
+            MagnitudeCube,      /// magCube(u)
+            MagnitudeSquareRoot /// magSqrt(u)
+        };
+
+
+    // - Member variables
+    
+        /// const reference to the Time object
+        const Time&         time_;
+
+        /// List to store all the point fields 
+        anyList             allFields_;
+
+        /// Map to store the last capture time of each field
+        wordMap<timeValue>  captureTime_;
+
+        /// 
+        bool                inSimulation_ = false;
+        
+        static const inline std::map<word, word> reservedFieldNames_
+        {
+            {"position", "realx3"},
+            {"one", "real"}
+        };
+        
+    // - Private member functions
+
+        uint32 pointFieldSize()
+        {
+            auto s = updatePoints();
+            return s.size();
+        }
+        
+        /// @brief Checks if a field needs to be updated.
+        /// @param name The name of the field to check.
+        /// @param forceUpdate Forces an update if true. Defaults to `false`.
+        /// @return `true` if the field needs updating or `forceUpdate` is true, `false` otherwise.
+        bool checkForUpdate(const word& name, bool forceUpdate = false);
+
+        /**
+         * @brief Updates and retrieves a point field of a specified type from the database.
+         * 
+         * This is a template function that updates and retrieves a point field of type `T`
+         * from the database. It checks if the field needs to be updated based on the last
+         * capture time or if a forced update is requested. If an update is necessary, it
+         * retrieves the latest data for the field.
+         * 
+         * @tparam T The type of the point field to update and retrieve. Must be a ValidFieldType.
+         * @param name The name of the field to update.
+         * @param forceUpdate A boolean flag indicating whether to force an update of the field
+         *                    regardless of its current state. Defaults to `false`.
+         * 
+         * @return A span of `T` representing the updated field data.
+         * 
+         * @throws message If the field cannot be found in the database or if there is a type mismatch.
+         */
+        template<ValidFieldType T>
+        span<T> updateField(const word& name, bool forceUpdate = false);
+        
+        /// @brief Creates a new real field or retrieves an existing one.
+        /// 
+        /// If a field with the given name already exists, it returns a span to that field.
+        /// If the field does not exist, it creates a new real field with the given name
+        /// and returns a span to the newly created field.
+        /// 
+        /// @param name The name of the field to create or retrieve.
+        /// @return span<real> of the field 
+        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.
+         */
+        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:
+    
+    // - constructors 
+
+        fieldsDataBase(const Time& time, bool inSimulation);
+
+        /// no copy constructor
+        fieldsDataBase(const fieldsDataBase&) = delete;
+        
+        /// no move constructor
+        fieldsDataBase(fieldsDataBase&&) = delete;
+        
+        /// no copy assignment
+        fieldsDataBase& operator=(const fieldsDataBase&) = delete;
+        
+        /// no move assignment
+        fieldsDataBase& operator=(fieldsDataBase&&) = delete;
+        
+        /// destructor
+        ~fieldsDataBase() = default;
+    
+    // - Public Access Functions
+        /// returns the current time
+        timeValue currentTime()const;
+
+        /// const ref to object Time
+        const Time& time()const
+        {
+            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(
+            const word& compoundName,
+            word& fieldName, 
+            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);
+
+        /**
+         * @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.
+         */
+        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
+         */
+        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);
+        
+        span<uint32> updateFieldUint32(
+            const word& name,
+            bool forceUpdate = false);
+    
+        /// Updates and retrieves a point field of any type from the database.
+        /// It returns types real, realx3 and realx4 only.
+        allPointFieldTypes updateFieldAll(
+            const word& compoundName, 
+            bool forceUpdate = false);
+
+        const pointStructure& pStruct()const;
+};
+
+} // nampespace pFlow
+
+#include "fieldsDataBaseTemplates.cpp"
+
+#endif //__fieldsDataBased_hpp__

src/PostProcessData/fieldsDataBase/fieldsDataBaseTemplates.cpp

@@ -0,0 +1,80 @@
+/*------------------------------- 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 __fieldsDataBaseTemplates_hpp__
+#define __fieldsDataBaseTemplates_hpp__
+
+#include "fieldsDataBase.hpp"
+
+template<pFlow::ValidFieldType T>
+inline
+pFlow::span<T> pFlow::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 constexpr( std::same_as<T, realx3>)
+            {
+                return updatePoints(forceUpdate);            
+            }
+            else
+            {
+                fatalErrorInFunction<< "Error in getting the type name of field: "<<
+                    name<<", with type name: "<< getTypeName<T>() <<endl;
+                fatalExit;
+                return span<T>(nullptr, 0);
+            }
+        }
+        else
+        {
+            const auto& pField = time_.lookupObject<pointField_D<T>>(name);
+            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
+            name, 
+            pField.activeValuesHost());
+        }
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<T>>(name);
+    
+    return span<T>(
+        field.data(), 
+        field.size());
+    
+}
+
+#endif //__fieldsDataBaseTemplates_hpp__

src/PostProcessData/method/GaussianDistribution.hpp

@@ -0,0 +1,107 @@
+/*------------------------------- 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 __GaussianDistribution_hpp__
+#define __GaussianDistribution_hpp__
+
+#include <vector>
+
+#include "typeInfo.hpp"
+#include "types.hpp"
+#include "span.hpp"
+
+
+namespace pFlow
+{
+
+class GaussianDistribution
+{
+private:
+
+    std::vector<real>      weight_;
+
+    real            variance_ = 1.0;
+
+    realx3          mean_ = {0,0,0};
+
+public:
+
+    TypeInfoNV("GaussianDistribution");
+
+
+    GaussianDistribution() = default;
+
+    GaussianDistribution(realx3 mean, real variance)
+    :
+        weight_(),
+        variance_(variance),
+        mean_(mean)
+    {}
+
+    GaussianDistribution(const GaussianDistribution&) = default;
+
+    GaussianDistribution(GaussianDistribution&&) = default;
+
+    GaussianDistribution& operator =(const GaussianDistribution&) = default;
+
+    GaussianDistribution& operator = (GaussianDistribution&&) = default;
+
+    ~GaussianDistribution()=default;
+
+    bool updateWeights(const realx3& center, const span<uint32>& indices,  const span<realx3>& points)
+    {
+        weight_.clear();
+        weight_.resize(indices.size());
+
+        real sum = 0.0;
+        for(uint32 i=0; i<indices.size(); i++)
+        {
+            auto x = points[indices[i]]-center;
+            auto f = exp(- dot(x,x)/(2*variance_));
+            weight_[i] = f;
+            sum += f;
+        }
+
+        for(auto& w: weight_)
+        {
+            w /= sum;
+        }   
+        return true;
+    }
+
+    bool updateWeights(uint32 n)
+    {
+        return false;
+    }
+
+    bool getWeight(uint32 i)const
+    {
+        return weight_[i];
+    }
+
+    span<real> getWeights()
+    {
+        return span<real>(weight_.data(), weight_.size());
+    }
+
+};
+} // namespace pFlow
+
+#endif // __GaussianDistribution_hpp__

src/PostProcessData/method/uniformDistribution.hpp

@@ -0,0 +1,86 @@
+/*------------------------------- 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 __uniformDistribution_hpp__
+#define __uniformDistribution_hpp__
+
+#include "types.hpp"
+#include "typeInfo.hpp"
+#include "span.hpp"
+
+namespace pFlow
+{
+
+class dictionary;
+
+class uniformDistribution
+{
+private:
+
+    std::vector<real>      weight_;
+    
+public:
+
+    // type info
+    TypeInfoNV("uniformDistribution");
+
+    uniformDistribution()
+    {}
+
+    uniformDistribution(const uniformDistribution&) = default;
+
+    uniformDistribution(uniformDistribution&&) = default;
+
+    uniformDistribution& operator=(const uniformDistribution&) = default;
+
+    uniformDistribution& operator=(uniformDistribution&&) = default;
+
+    ~uniformDistribution()=default;
+
+
+    bool updateWeights(const realx3& center, const span<realx3>& points)
+    {
+        uint32 n = max(points.size(), 1u);
+        weight_.assign(n, 1.0/n);
+        return true;
+    }
+
+    bool updateWeights(uint32 n)
+    {
+        n = max(n, 1u);
+        weight_.assign(n, 1.0/n);
+        return true;
+    }
+
+    real getWeight(uint32 i)const
+    {
+        return weight_[i];
+    }
+
+    span<real> getWeights()
+    {
+        return span<real>(weight_.data(), weight_.size());
+    }
+
+};
+
+}
+
+#endif //__uniformDistribution_hpp__

src/PostProcessData/operation/PostprocessOperationSum.cpp

@@ -0,0 +1,69 @@
+#include "PostprocessOperationSum.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+#include "fieldFunctions.hpp"
+
+pFlow::PostprocessOperationSum::PostprocessOperationSum
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(opDict, regPoints, fieldsDB)
+{
+    if( fieldType() == getTypeName<real>() )
+    {
+        processedRegField_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, real(0)));
+    }
+    else if( fieldType() == getTypeName<realx3>() )
+    {
+        processedRegField_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx3(0)));
+    }
+    else if( fieldType() == getTypeName<realx4>() )
+    {
+        processedRegField_ = makeUnique<processedRegFieldType>(
+            regionField(processedFieldName(), regPoints, realx4(0)));
+    }
+    else
+    {
+        fatalErrorInFunction<<" in dictionary "<< opDict.globalName()
+            << " field type is not supported for sum operation"
+            << " field type is "<< fieldType()
+            << endl;
+        fatalExit;
+    }
+}
+
+bool pFlow::PostprocessOperationSum::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();
+
+
+    std::visit([&](auto&& field)->processedRegFieldType
+        { 
+            return executeSumOperation(
+            procName,
+            field,
+            regP,
+            dbVol,
+            weights,
+            phi,
+            mask);
+        },
+        allField);
+   
+    return true;
+}

src/PostProcessData/operation/PostprocessOperationSum.hpp

@@ -0,0 +1,76 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __PostprocessOperationSum_hpp__
+#define __PostprocessOperationSum_hpp__
+
+#include <variant>
+#include <vector>
+
+#include "postprocessOperation.hpp"
+#include "regionField.hpp"
+#include "includeMask.hpp"
+
+namespace pFlow
+{
+
+
+class PostprocessOperationSum
+:
+    public postprocessOperation
+{
+private:
+       
+    using processedRegFieldType = std::variant
+    <
+        regionField<real>,
+        regionField<realx3>,
+        regionField<realx4>
+    >;
+    
+    /// Pointer to the include mask used for masking operations.
+    uniquePtr<processedRegFieldType>  processedRegField_ = nullptr;
+
+public:
+
+    TypeInfo("PostprocessOperation<sum>");
+
+    PostprocessOperationSum(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    ~PostprocessOperationSum() override = default;
+
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationSum,
+        dictionary
+    );
+
+    bool execute(const std::vector<span<real>>& weights) override;
+
+}; 
+
+
+}
+
+#endif //__PostprocessOperation_hpp__

src/PostProcessData/operation/fieldFunctions.hpp

@@ -0,0 +1,117 @@
+/*------------------------------- 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 __fieldFunctions_hpp__
+#define __fieldFunctions_hpp__
+
+#include <vector>
+
+#include "span.hpp"
+#include "regionPoints.hpp"
+#include "regionField.hpp"
+#include "includeMask.hpp"
+
+namespace pFlow
+{
+
+template<typename T>
+regionField<T> executeSumOperation
+(
+    const word&                     regFieldName,
+    const span<T>&                  field, 
+    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{});
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto vols = regPoints.volumes();
+        auto w = weights[reg];
+        T sum{};
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( mask( index ))
+            {
+                sum += w[n] * field[index]* phi[index];
+            }
+            n++;
+        }
+        if(devideByVol)
+        {
+            processedField[reg] = sum/vols[reg];
+        }
+        else
+        {
+            processedField[reg] = sum;
+        }
+        
+    }
+
+    return processedField;
+}
+
+template<typename T>
+regionField<T> executeAverageOperation
+(
+    const word&                    regFieldName,
+    const span<T>&                 field, 
+    const regionPoints&            regPoints, 
+    const std::vector<span<real>>& weights,
+    const span<real>&              phi,
+    const includeMask::Mask&       mask
+)
+{
+    regionField<T> processedField(regFieldName, regPoints, T{});
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto w = weights[reg];
+        T sumNum{};
+        real sumDen{};
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( mask( index ))
+            {
+                sumNum += w[n] * field[index]* phi[index];
+            }
+            sumDen += w[n] * phi[index];
+            n++;
+        }
+
+        sumDen = max(sumDen, smallValue);
+        processedField[reg] = sumNum/sumDen;
+        
+    }
+
+    return processedField;
+}
+
+} // namespace pFlow
+
+#endif //__fieldFunctions_hpp__

src/PostProcessData/operation/includeMask/IncludeMask.hpp

@@ -0,0 +1,192 @@
+/*------------------------------- 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 __IncludeMask_hpp__
+#define __IncludeMask_hpp__
+
+#include "dictionary.hpp"
+#include "includeMask.hpp"
+#include "fieldsDataBase.hpp"
+#include "maskOperations.hpp"
+#include "Time.hpp"
+
+
+namespace pFlow
+{
+
+
+template<typename T, typename Operator>
+class IncludeMask
+:
+	public includeMask
+{
+public:
+
+    using Mask = typename includeMask::Mask;
+
+private:
+	
+    std::vector<bool>       mask_;
+
+    Operator                operator_;
+
+    word                    fieldName_;
+    
+    timeValue               lastUpdated_ = -1;
+
+    bool updateMask()
+    {
+        timeValue t = database().currentTime();
+        
+        if( equal( t, lastUpdated_)) return true;
+        
+        span<T> s;
+        if constexpr (std::is_same_v<T,real>)
+        {
+            s = database().updateFieldReal(fieldName_);
+        }
+        else if constexpr ( std::is_same_v<T,realx3>)
+        {
+            s = database().updateFieldRealx3(fieldName_);
+        }
+        else if constexpr( std::is_same_v<T,realx4>)
+        {
+            s = database().updateFieldRealx4(fieldName_);   
+        }
+        else
+        {
+            fatalErrorInFunction<<"Type "<< getTypeName<T>()
+                <<" is not supported for  IncludeMask for field "
+                << fieldName_ << endl;
+            return false;
+        }
+        
+        
+        mask_.resize(s.size());
+
+        for(uint32 i=0; i<s.size(); i++)
+        {
+            mask_[i] = operator_(s[i]);
+        }
+
+        lastUpdated_ = t ;
+        
+        return true;
+    }
+
+    static
+    word operatorName()
+    {
+        return Operator::TYPENAME();
+    }
+
+public:
+
+	TypeInfoTemplate12("IncludeMask", T, Operator);
+
+	IncludeMask(
+		const dictionary& dict, 
+        fieldsDataBase& feildsDB)
+	:
+		includeMask(dict, feildsDB),
+		operator_(dict.subDict(operatorName()+"Info")),
+        fieldName_(
+            dict.subDict
+            (
+                operatorName()+"Info"
+            ).getVal<word>("field"))
+	{}
+
+	add_vCtor(
+		includeMask,
+		IncludeMask,
+		dictionary);
+
+    Mask getMask() override
+    {
+        updateMask();
+        return Mask(mask_);
+    }
+	
+};
+
+
+template<typename T>
+class IncludeMask<T,allOp<T>>
+:
+	public includeMask
+{
+public:
+
+    using Mask = typename includeMask::Mask;
+
+private:
+	
+    std::vector<bool>       mask_;
+
+    timeValue               lastUpdated_ = -1;
+
+    bool updateMask()
+    {
+        timeValue t = database().currentTime();
+        
+        if( equal( t, lastUpdated_)) return true;
+        
+        span<realx3>  s = database().updatePoints();        
+        mask_.resize(s.size(), true);
+
+        lastUpdated_ = t ;
+        
+        return true;
+    }   
+
+public:
+
+	TypeInfoTemplate12("IncludeMask", T, allOp<int8>);
+
+	IncludeMask(
+		const dictionary& opDict, 
+        fieldsDataBase& feildsDB)
+	:
+		includeMask(opDict, feildsDB)
+	{
+        span<realx3>  s = database().updatePoints();        
+        mask_.resize(s.size(), true);
+	}
+
+	add_vCtor(
+		includeMask,
+		IncludeMask,
+		dictionary);
+
+    Mask getMask()override
+    {
+        updateMask();
+        return Mask(mask_);
+    }
+	
+};
+
+
+} // pFlow
+
+#endif //__IncludeMask_hpp__
+
+

src/PostProcessData/operation/includeMask/IncludeMasks.cpp

@@ -0,0 +1,50 @@
+/*------------------------------- 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 "IncludeMask.hpp"
+
+// real
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareOne<pFlow::real, pFlow::lessThanOp> >;
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareOne<pFlow::real, pFlow::lessThanEqOp> >;
+
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareOne<pFlow::real, pFlow::greaterThanOp> >;
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareOne<pFlow::real, pFlow::greaterThanEqOp> >;
+
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareOne<pFlow::real, pFlow::equalOp> >;
+
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareTwo<pFlow::real, pFlow::betweenOp> >;
+template class pFlow::IncludeMask<pFlow::real, pFlow::compareTwo<pFlow::real, pFlow::betweenEqOp> >;
+
+template class pFlow::IncludeMask<pFlow::real, pFlow::allOp<pFlow::real>>;
+
+// realx3
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareOne<pFlow::realx3, pFlow::lessThanOp> >;
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareOne<pFlow::realx3, pFlow::lessThanEqOp> >;
+
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareOne<pFlow::realx3, pFlow::greaterThanOp> >;
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareOne<pFlow::realx3, pFlow::greaterThanEqOp> >;
+
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareOne<pFlow::realx3, pFlow::equalOp> >;
+
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareTwo<pFlow::realx3, pFlow::betweenOp> >;
+template class pFlow::IncludeMask<pFlow::realx3, pFlow::compareTwo<pFlow::realx3, pFlow::betweenEqOp> >;
+
+
+// realx4

src/PostProcessData/operation/includeMask/includeMask.cpp

@@ -0,0 +1,91 @@
+/*------------------------------- 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 "includeMask.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+
+
+pFlow::includeMask::includeMask
+(
+	const dictionary& dict, 
+	fieldsDataBase& fieldDB
+)
+:
+	database_(fieldDB)
+{}
+
+
+pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
+(
+    const dictionary& opDict, 
+    fieldsDataBase& feildsDB
+)
+{
+	word mask = opDict.getValOrSet<word>("includeMask", "all");
+	word fieldType;
+	if( mask != "all")
+	{
+		auto& maskDict = opDict.subDict(mask+"Info");
+		word maskField = maskDict.getVal<word>("field");
+		
+		if( !feildsDB.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, mask);
+	
+	if( dictionaryvCtorSelector_.search(method) )
+	{
+		auto objPtr = 
+			dictionaryvCtorSelector_[method]
+			(opDict, feildsDB);
+		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

@@ -0,0 +1,115 @@
+/*------------------------------- 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 __includeMask_hpp__
+#define __includeMask_hpp__
+
+#include <vector> 
+
+#include "virtualConstructor.hpp"
+
+namespace pFlow
+{
+
+
+class fieldsDataBase;
+class dictionary;
+
+class includeMask
+{
+public:
+
+	class Mask
+	{
+		const std::vector<bool>& mask_;
+
+		public:
+
+		Mask(const std::vector<bool>& msk)
+		:
+			mask_(msk)
+		{}
+
+		Mask(const Mask&) = default;
+
+		Mask(Mask&&) = default;
+				
+		~Mask()=default;
+
+		auto size()const
+		{
+			return mask_.size();
+		}
+
+		bool operator()(uint32 i)const
+		{
+			return mask_[i];
+		}
+	};
+
+private:
+	
+    fieldsDataBase&     database_;
+	
+public:
+
+	TypeInfo("includeMask");
+
+	includeMask(const dictionary& opDict, fieldsDataBase& feildsDB);
+
+	virtual ~includeMask() = default;
+
+    create_vCtor
+    (
+        includeMask,
+        dictionary,
+        (
+            const dictionary& opDict, fieldsDataBase& feildsDB
+        ),
+        (opDict, feildsDB)		
+    );
+	 
+	const fieldsDataBase& database()const
+	{
+		return database_;
+	}
+
+    fieldsDataBase& database()
+    {
+        return database_;
+    }
+
+    virtual 
+    Mask getMask()= 0;
+
+	static 
+	uniquePtr<includeMask> create(
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB);
+	
+};
+
+
+
+} // pFlow
+
+#endif //__IncludeMask_hpp__
+
+

src/PostProcessData/operation/includeMask/maskOperations.hpp

@@ -0,0 +1,162 @@
+
+#ifndef __maskOperation_hpp__
+#define __maskOperation_hpp__
+
+#include "types.hpp"
+#include "dictionary.hpp"
+
+namespace pFlow
+{
+
+template<typename T>
+struct greaterThanOp
+{
+	TypeInfoNV("greaterThan");
+
+ 	inline
+ 	bool operator()(const T &compVal, const T &val) const {
+   		return val > compVal; }
+};
+
+template<typename T>
+struct greaterThanEqOp
+{
+	TypeInfoNV("greaterThanEq");
+
+ 	inline
+ 	bool operator()(const T &compVal, const T &val) const {
+   		return val >= compVal; }
+};
+
+template<typename T>
+struct lessThanOp
+{
+	TypeInfoNV("lessThan");
+
+	inline
+	bool operator()(const T &compVal, const T &val) const {
+		return val < compVal; }
+};
+
+template<typename T>
+struct lessThanEqOp
+{
+	TypeInfoNV("lessThanEq");
+
+	inline
+	bool operator()(const T &compVal, const T &val) const {
+		return val <= compVal; }
+};
+
+template<typename T>
+struct equalOp
+{
+	TypeInfoNV("equal");
+
+	inline
+	bool operator()(const T &compVal, const T &val) const {
+		return equal(val , compVal); }
+};
+
+
+template<typename T>
+struct betweenOp
+{
+	TypeInfoNV("between");
+
+	inline
+	bool operator()(const T &compVal1, const T &compVal2 ,const T &val) const {
+		return val>compVal1 && val<compVal2; }
+};
+
+
+template<typename T>
+struct betweenEqOp
+{
+	TypeInfoNV("betweenEq");
+
+	inline
+	bool operator()(const T &compVal1, const T &compVal2 ,const T &val) const {
+		return val>=compVal1 && val<=compVal2; }
+};
+
+template<typename T>
+struct allOp
+{
+	TypeInfoNV("all");
+
+	inline
+	bool operator()() const {return true; }
+};
+
+
+
+template<typename T, template<class> class Operator>
+class compareOne
+{
+public:
+
+	using opertorType = Operator<T>;
+
+protected:
+	T 				compValue_{};
+	opertorType 	operator_{};
+public:
+
+	TypeInfoNV(Operator<T>::TYPENAME());
+
+	compareOne(const dictionary& dict)
+	:
+		compValue_(dict.getVal<T>("value"))
+	{}
+
+	bool operator()(const T& value)const
+	{
+		return operator_(compValue_, value);
+	}
+};
+
+template<typename T, template<class> class Operator>
+class compareTwo
+{
+public:
+	using opertorType = Operator<T>;
+protected:
+	T 				compValue1_;
+	T  				compValue2_;
+	opertorType 	operator_{};
+public:
+	
+	TypeInfoNV(opertorType::TYPENAME());
+
+	compareTwo(const dictionary& dict)
+	:
+		compValue1_(dict.getVal<T>("value1")),
+		compValue2_(dict.getVal<T>("value2"))
+	{}
+
+	bool operator()(const T& value)const
+	{
+		return operator_(compValue1_, compValue2_, value);
+	}
+};
+
+template<typename T, typename Operator>
+class compareZero
+{
+protected:
+	Operator 	operator_{};
+public:
+
+	TypeInfoNV(Operator::TYPENAME());
+	compareZero(const dictionary& dict);
+
+	bool operator()(const T& value) const
+	{
+		return operator_();
+	}
+};
+
+}
+
+#endif //__maskOperation_hpp__

src/PostProcessData/operation/postprocessOperation.cpp

@@ -0,0 +1,101 @@
+/*------------------------------- 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 "postprocessOperation.hpp"
+#include "regionPoints.hpp"
+#include "fieldsDataBase.hpp"
+
+pFlow::postprocessOperation::postprocessOperation
+(
+    const dictionary &opDict, 
+    const regionPoints& regPoints,
+    fieldsDataBase &fieldsDB
+)
+:
+    operationDict_(opDict),
+    threshold_
+    (
+        opDict.getValOrSet<int>("threshold", 1)
+    ),
+    divideByVolume_
+    (
+        opDict.getValOrSet<Logical>("dividedByVolume", Logical(false))
+    ),
+    regionPoints_
+    (
+        regPoints
+    ),
+    database_
+    (
+        fieldsDB
+    ),
+    fieldName_
+    (
+        opDict.getValOrSet<word>("field", "one")
+    ),
+    phiFieldName_
+    (
+        opDict.getValOrSet<word>("phi", "one")
+    ),
+    includeMask_
+    (
+        includeMask::create(opDict, fieldsDB)
+    )
+{
+
+    if(!fieldsDB.getPointFieldType(fieldName_, fieldType_))
+    {
+        fatalErrorInFunction;
+        fatalExit;
+    }
+}
+
+pFlow::uniquePtr<pFlow::postprocessOperation> 
+    pFlow::postprocessOperation::create
+    (
+        const dictionary &opDict, 
+        const regionPoints& regPoints,
+        fieldsDataBase &fieldsDB
+    )
+{
+    word func = opDict.getVal<word>("function");
+    word method = angleBracketsNames("PostprocessOperation", func);
+    if( dictionaryvCtorSelector_.search(method) )
+    {
+        REPORT(3)<<"Operation "<< Green_Text(opDict.name())<<" with function "<< Green_Text(func)<<endl;
+        auto objPtr = 
+            dictionaryvCtorSelector_[method]
+            (opDict, regPoints, fieldsDB);
+        return objPtr;
+    }
+    else
+    {
+        printKeys
+        ( 
+            fatalError << "Ctor Selector "<< 
+            method << " dose not exist. \n"
+            <<"Avaiable ones are: \n\n"
+            ,
+            dictionaryvCtorSelector_
+        );
+        fatalExit;
+        return nullptr;
+    }
+}

src/PostProcessData/operation/postprocessOperation.hpp

@@ -0,0 +1,167 @@
+/*------------------------------- 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 __postprocessOperation_hpp__
+#define __postprocessOperation_hpp__
+
+#include "virtualConstructor.hpp"
+#include "Logical.hpp"
+#include "dictionary.hpp"
+#include "span.hpp"
+#include "includeMask.hpp"
+
+
+namespace pFlow
+{
+
+class fieldsDataBase;
+class regionPoints;
+
+/*!
+ *  @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:
+
+    using Mask = typename includeMask::Mask;
+
+private:
+
+    /// Dictionary containing operation-specific parameters.
+    dictionary              operationDict_;
+
+    /// This Threshold is used to exclude the regions which contain
+    /// fewer than this value. 
+    uint32                  threshold_;
+    
+    /// Logical flag to determine if the result is divided by volume.
+    Logical                 divideByVolume_;
+
+    /// Reference to the region points used in the operation.
+    const regionPoints&     regionPoints_;
+
+    /// Reference to the fields database containing field data.
+    fieldsDataBase&         database_;
+
+    /// Name of the field to be processed.
+    word                    fieldName_;
+
+    /// Type of the field to be processed.
+    word                    fieldType_;
+
+    /// Name of the phi field to be processed.
+    word                     phiFieldName_;
+
+    /// Pointer to the include mask used for masking operations.
+    uniquePtr<includeMask>  includeMask_ = nullptr;
+    
+public:
+
+
+    TypeInfo("postprocessOperation");
+
+    postprocessOperation(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB );
+
+    virtual ~postprocessOperation()=default;
+
+    create_vCtor(
+        postprocessOperation,
+        dictionary,
+        (
+            const dictionary& opDict,
+            const regionPoints& regPoints,
+            fieldsDataBase& fieldsDB
+        ),
+        (opDict, regPoints, fieldsDB));
+    
+    const regionPoints& regPoints()const
+    {
+        return regionPoints_;
+    }
+
+    const fieldsDataBase& database()const
+    {
+        return database_;
+    }
+
+    fieldsDataBase& database()
+    {
+        return database_;
+    }
+
+    word processedFieldName()const
+    {
+        return operationDict_.name();
+    }
+    
+    const word& fieldName()const
+    {
+        return fieldName_;
+    }
+
+    const word& fieldType()const
+    {
+        return fieldType_;
+    }
+
+    const word& phiFieldName()const
+    {
+        return phiFieldName_;
+    }
+
+    const dictionary& operationDict()const
+    {
+        return operationDict_;
+    }
+    const uint32 threshold()const
+    {
+        return threshold_;
+    }
+
+    bool divideByVolume()const
+    {
+        return divideByVolume_();
+    }
+
+    Mask getMask()
+    {
+        return includeMask_().getMask();
+    }
+
+
+    virtual 
+    bool execute(const std::vector<span<real>>& weights) = 0;
+
+    static
+    uniquePtr<postprocessOperation> create(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+};
+
+}
+
+#endif //__postprocessOperation_hpp__

src/PostProcessData/postprocessComponent/PostprocessComponent.cpp

@@ -0,0 +1,109 @@
+#include "PostprocessComponent.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.
+
+-----------------------------------------------------------------------------*/
+
+template<typename RegionType, typename ProcessMethodType>
+pFlow::PostprocessComponent<RegionType,ProcessMethodType>::PostprocessComponent
+(
+    const dictionary& dict,
+    fieldsDataBase& fieldsDB,
+    const baseTimeControl& defaultTimeControl
+)
+:
+    postprocessComponent(dict, fieldsDB, defaultTimeControl),
+    regionPointsPtr_
+    (
+        makeUnique<RegionType>(dict, fieldsDB)
+    ),
+    regionsProcessMethod_
+    (
+        regionPointsPtr_().size()
+    ),
+    operationDicts_(readDictList("operations", dict))
+{
+
+    for(auto& opDict:operationDicts_)
+    {
+        operatios_.push_back
+        (
+            postprocessOperation::create
+            (
+                opDict,
+                regionPointsPtr_(), 
+                this->database()
+            )
+        );
+    }
+}
+
+
+template <typename RegionType, typename ProcessMethodType>
+bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::execute
+(
+    const timeInfo &ti, 
+    bool forceUpdate
+)
+{
+    
+    if( !timeControl().eventTime(ti))
+    {
+        executed_ = false;
+        return true;
+    }
+
+    // update processing methods
+    
+    auto centers = this->regPoints().centers();
+    const uint32 n     = centers.size();
+    auto points  = this->database().updatePoints();
+    
+    for(uint32 i=0; i<n; i++)
+    {
+        auto indices = this->regPoints().indices(i);
+        regionsProcessMethod_[i].updateWeights(
+            centers[i],
+            indices,
+            points);
+    }
+
+    std::vector<span<real>> weights(n);
+
+    for(uint32 i=0; i<n; i++)
+    {
+        // get the span of weight of each region 
+        weights[i] = regionsProcessMethod_[i].getWeights();
+    }
+    
+    for(auto& op:operatios_)
+    {
+        if( op->execute(weights) )
+        {
+            fatalErrorInFunction
+                <<"error occured in executing operatoin defined in dict "
+                << op->operationDict()
+                <<endl;
+            return false;
+        }
+    }
+
+    executed_ = true;
+
+    return true;
+}

src/PostProcessData/postprocessComponent/PostprocessComponent.hpp

@@ -0,0 +1,117 @@
+/*------------------------------- 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 __PostProcessComponent_hpp__
+#define __PostProcessComponent_hpp__
+
+#include <vector>
+
+#include "ListPtr.hpp"
+#include "List.hpp"
+#include "postprocessComponent.hpp"
+#include "postprocessOperation.hpp"
+#include "dictionaryList.hpp"
+#include "fieldsDataBase.hpp"
+#include "regionPoints.hpp"
+#include "regionField.hpp"
+
+namespace pFlow
+{
+
+template<typename RegionType, typename ProcessMethodType>
+class PostprocessComponent
+:
+    public postprocessComponent
+{
+private:
+
+
+    ListPtr<postprocessOperation> operatios_; 
+
+    /// Region type for selecting a subset of particles for processing
+    uniquePtr<RegionType>         regionPointsPtr_;
+    
+    /// Method for processing the selected particles data 
+    std::vector<ProcessMethodType>  regionsProcessMethod_;
+
+    bool    executed_{false};
+
+    dictionaryList                  operationDicts_;
+
+protected:
+    
+    std::vector<ProcessMethodType>& regionProecessMethod()
+    {
+        return regionsProcessMethod_;
+    }
+
+public:
+
+    // type info 
+    TypeInfoTemplate12(
+        "PostprocessComponent",
+        RegionType,
+        ProcessMethodType);
+
+    PostprocessComponent(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl);
+    
+    ~PostprocessComponent() override = default;
+
+    // add the virtual constructor
+    add_vCtor(
+        postprocessComponent,
+        PostprocessComponent,
+        dictionary
+    );   
+
+    word name()const override
+    {
+        return operationDicts_.parrent().name();
+    }
+
+    regionPoints& regPoints() override
+    {
+        return static_cast<regionPoints&>(regionPointsPtr_());
+    }
+
+    const regionPoints& regPoints()const override
+    {
+        return static_cast<const regionPoints&>(regionPointsPtr_());
+    }
+
+    bool execute(const timeInfo& ti, bool forceUpdate = false) override;
+
+    bool executed()const override
+    {
+        return executed_;
+    }
+
+
+
+};
+
+}
+
+#include "PostprocessComponent.cpp"
+
+#endif

src/PostProcessData/postprocessComponent/PostprocessComponentGaussian.hpp

@@ -0,0 +1,73 @@
+/*------------------------------- 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 __PostprocessComponentGaussian_hpp__
+#define __PostprocessComponentGaussian_hpp__
+
+#include "PostprocessComponent.hpp"
+#include "GaussianDistribution.hpp"
+
+namespace pFlow
+{
+
+template<typename RegionType>
+class PostprocessComponentGaussian
+:
+    public PostprocessComponent<RegionType, GaussianDistribution>
+{
+public:
+
+    /// type info
+    TypeInfoTemplate12("PostprocessComponent", RegionType, GaussianDistribution);
+
+    PostprocessComponentGaussian
+    (
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl
+    )
+    :
+        PostprocessComponent<RegionType,GaussianDistribution>(dict, fieldsDB, defaultTimeControl)
+    {
+        /// initializes the Gaussian distribution for all elements of region
+        //const uint32 n = this->regPoints().size();
+        auto d = this->regPoints().eqDiameters();
+        auto c = this->regPoints().centers();
+        auto& regs = this->regionProecessMethod();
+        const uint32 n = d.size();
+        for(uint32 i=0; i<n; i++)
+        {
+            regs[i] = GaussianDistribution(c[i], pow(d[i],2));
+        }
+    }
+
+    // add the virtual constructor
+    add_vCtor(
+        postprocessComponent,
+        PostprocessComponentGaussian,
+        dictionary
+    );
+
+};
+
+}
+
+
+#endif //__PostprocessComponentGaussian_hpp__

src/PostProcessData/postprocessComponent/PostprocessComponents.cpp

@@ -0,0 +1,34 @@
+/*------------------------------- 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 "PostprocessComponentGaussian.hpp"
+
+// region types 
+#include "sphereRegionPoints.hpp"
+
+
+// methods 
+#include "GaussianDistribution.hpp"
+#include "uniformDistribution.hpp"
+
+
+template class pFlow::PostprocessComponentGaussian<pFlow::sphereRegionPoints>;
+
+

src/PostProcessData/postprocessComponent/particleProbePostprocessComponent.cpp

@@ -0,0 +1,26 @@
+#include "particleProbePostprocessComponent.hpp"
+
+pFlow::particleProbePostprocessComponent::particleProbePostprocessComponent
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDB, 
+    const baseTimeControl &defaultTimeControl
+)
+:
+    postprocessComponent(dict, fieldsDB, defaultTimeControl),
+    regionPointsPtr_
+    (
+        makeUnique<centerPointsRegionPoints>(dict, fieldsDB)
+    ),
+    name_(dict.name())
+{}
+
+bool pFlow::particleProbePostprocessComponent::execute
+(
+    const timeInfo &ti, 
+    bool forceExecute
+)
+{
+    
+    return false;
+}

src/PostProcessData/postprocessComponent/particleProbePostprocessComponent.hpp

@@ -0,0 +1,94 @@
+/*------------------------------- 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 __particleProbePostprocessComponent_hpp__
+#define __particleProbePostprocessComponent_hpp__
+
+#include "postprocessComponent.hpp"
+#include "fieldsDataBase.hpp"
+#include "centerPointsRegionPoints.hpp"
+#include "regionField.hpp"
+
+namespace pFlow
+{
+
+class particleProbePostprocessComponent
+:
+    public postprocessComponent
+{
+private:
+
+    using processedRegFieldType = std::variant
+    <
+        regionField<real>,
+        regionField<realx3>,
+        regionField<realx4>
+    >;
+
+    bool                                executed_{false};
+
+    uniquePtr<centerPointsRegionPoints> regionPointsPtr_;
+
+    uniquePtr<processedRegFieldType>    processedField_ = nullptr;
+
+    word                                name_;
+
+public:
+
+    TypeInfo("PostprocessComponent<centerPoints,particleProbe>");
+
+    particleProbePostprocessComponent
+    (
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl
+    );
+    
+    ~particleProbePostprocessComponent()override = default;
+
+     
+    word name()const override
+    {
+        return name_;
+    }
+    
+    regionPoints& regPoints() override
+    {
+        return regionPointsPtr_();
+    }
+
+     
+    const regionPoints& regPoints() const override
+    {
+        return regionPointsPtr_();
+    }
+
+    bool execute(const timeInfo& ti, bool forceExecute = false) override;
+
+    bool executed() const override
+    {
+        return executed_;
+    }
+
+};
+    
+}
+
+#endif //__particleProbePostprocessComponent_hpp__

src/PostProcessData/postprocessComponent/postprocessComponent.cpp

@@ -0,0 +1,55 @@
+#include "postprocessComponent.hpp"
+#include "fieldsDataBase.hpp"
+#include "Time.hpp"
+
+
+pFlow::postprocessComponent::postprocessComponent
+( 
+    const dictionary &dict,
+    fieldsDataBase &fieldsDB, 
+    const baseTimeControl &defaultTimeControl
+)
+:
+    timeControl_(
+        dict,
+        defaultTimeControl,
+        baseTimeControl(
+            fieldsDB.time().startTime(),
+            fieldsDB.time().endTime(),
+            fieldsDB.time().saveInterval())
+    ),
+    fieldsDataBase_(fieldsDB)
+{
+}
+
+
+pFlow::uniquePtr<pFlow::postprocessComponent> pFlow::postprocessComponent::create
+(
+    const dictionary& dict,
+    fieldsDataBase& fieldsDB,
+    const baseTimeControl& defaultTimeControl
+)
+{
+    word method = dict.getVal<word>("processMethod");
+    word region = dict.getVal<word>("processRegion");
+    
+    auto compModel = angleBracketsNames2("PostprocessComponent", region, method);
+
+    if( dictionaryvCtorSelector_.search(compModel) )
+    {
+        REPORT(2)<<"Creating postprocess component "<< Green_Text(compModel)<<" ..."<<endl;
+        return dictionaryvCtorSelector_[compModel](dict, fieldsDB, defaultTimeControl);
+    }
+    else
+    {
+        printKeys
+		( 
+			fatalError << "Ctor Selector "<< Yellow_Text(compModel) << " dose not exist. \n"
+			<<"Avaiable ones are: \n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+        return nullptr;
+    }
+}

src/PostProcessData/postprocessComponent/postprocessComponent.hpp

@@ -0,0 +1,114 @@
+/*------------------------------- 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 __postprocessComponent_hpp__
+#define __postprocessComponent_hpp__
+
+#include "postprocessTimeControl.hpp"
+#include "dictionary.hpp"
+#include "virtualConstructor.hpp"
+
+namespace pFlow
+{
+
+class fieldsDataBase;
+class regionPoints;
+class dictionary;
+
+class postprocessComponent
+{
+private:
+
+    postprocessTimeControl  timeControl_;
+
+    fieldsDataBase&         fieldsDataBase_;
+
+public:
+
+    /// type info
+    TypeInfo("postprocessComponent");
+
+    postprocessComponent( 
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl);
+    
+    virtual ~postprocessComponent() = default;
+
+
+    create_vCtor
+    (
+        postprocessComponent,
+        dictionary,
+        (
+            const dictionary& dict, 
+            fieldsDataBase& fieldsDB, 
+            const baseTimeControl& defaultTimeControl
+        ),
+        (dict, fieldsDB, defaultTimeControl)
+    );
+
+    /// @brief  return the time control for this post-process component
+    /// @return reference to the time control object
+    const postprocessTimeControl& timeControl()const
+    {
+        return timeControl_;
+    }
+    
+    /// @brief  return the fields database
+    /// @return const reference to the fields database
+    const fieldsDataBase& database()const
+    {
+        return fieldsDataBase_;
+    }
+
+    /// @brief  return the fields database
+    /// @return reference to the fields database
+    fieldsDataBase& database()
+    {
+        return fieldsDataBase_;
+    }
+
+    virtual 
+    word name()const = 0;
+
+    virtual 
+    regionPoints& regPoints() = 0;
+
+    virtual 
+    const regionPoints& regPoints() const = 0;
+
+    virtual 
+    bool execute(const timeInfo& ti, bool forceExecute = false) = 0;
+
+    virtual 
+    bool executed()const = 0;
+
+    static
+    uniquePtr<postprocessComponent> create(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl);
+
+};
+
+} // namespace pFlow
+
+#endif // __postprocessComponent_hpp__

src/PostProcessData/postprocessData.cpp

@@ -0,0 +1,99 @@
+
+/*------------------------------- 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 "List.hpp"
+#include "systemControl.hpp"
+#include "postprocessData.hpp"
+#include "fileDictionary.hpp"
+#include "postprocessComponent.hpp"
+
+pFlow::postprocessData::postprocessData(const systemControl &control)
+:
+    auxFunctions(control),
+    time_(control.time()),
+    fieldsDataBase_(control.time(), true),
+    dict_
+    (
+        objectFile
+        (
+            "postprocessDataDict",
+            control.settings().path(),
+            objectFile::READ_IF_PRESENT,
+            objectFile::WRITE_NEVER
+        )
+    ),
+    componentsDicts_(readDictList("components", dict_))
+{
+    // if dictionary is not provided, no extra action is required. 
+    if( !dict_.fileExist() )
+    {
+        return;
+    }
+
+    activeInSimulation_ = dict_.getValOrSet<Logical>(
+        "activeInSimulation", 
+        Logical{true});
+
+    if(dict_.containsDictionay("defaultTimeControl"))
+    {
+        defaultTimeControlPtr_ = 
+            makeUnique<baseTimeControl>(
+                dict_.subDict("defaultTimeControl"),
+                "execution");
+    }
+    else
+    {
+        // default time control from settings
+        defaultTimeControlPtr_ = makeUnique<baseTimeControl>(
+            control.time().startTime(), 
+            control.time().endTime(),
+            control.time().saveInterval(),
+            "execution");
+    }
+    
+    for(auto& compDict:componentsDicts_)
+    {
+        postprocesses_.push_back( postprocessComponent::create(
+            compDict, 
+            fieldsDataBase_, 
+            defaultTimeControlPtr_() ));
+    }
+
+}
+
+bool pFlow::postprocessData::execute() 
+{
+    
+    const auto& ti = time_.TimeInfo();
+
+    for(auto& component:postprocesses_)
+    {
+        if(!component->execute(ti))
+        {
+            fatalErrorInFunction
+                <<"Error occured in executing postprocess component: "
+                <<component->name()<<endl;
+            return false;
+        }
+    }
+
+    return true;
+}

src/PostProcessData/postprocessData.hpp

@@ -0,0 +1,103 @@
+/*------------------------------- 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 __postprocessData_hpp__
+#define __postprocessData_hpp__
+
+#include "auxFunctions.hpp"
+#include "Logical.hpp"
+#include "ListPtr.hpp"
+#include "fileDictionary.hpp"
+#include "baseTimeControl.hpp"
+#include "fieldsDataBase.hpp"
+#include "postprocessComponent.hpp"
+#include "dictionaryList.hpp"
+
+namespace pFlow 
+{
+
+class systemControl;
+class Time;
+class timeInfo;
+
+/**
+ * @class postprocessData
+ * @brief An interface class for handling post-processing of simulation data.
+ * 
+ * This class provides methods and utilities to process and analyze
+ * simulation data during/after simulation.
+ */
+class postprocessData
+:
+  public auxFunctions
+{
+    /// Indicates if a post-processing is active during simulatoin 
+    Logical                 activeInSimulation_{false};
+
+    /// a list of active post-process components 
+    ListPtr<postprocessComponent>    postprocesses_;
+
+    /// const ref to Time
+    const Time&           time_;
+
+    /// Database for all the points fields on the host 
+    fieldsDataBase        fieldsDataBase_;
+
+    /// file dictionary that is constructed from the file (postProcessDataDict)
+    fileDictionary        dict_;
+
+    /// list of dictionaries for postprocess components 
+    dictionaryList        componentsDicts_;
+
+    /// @brief  default time control that can be used for all post-process components
+    uniquePtr<baseTimeControl> defaultTimeControlPtr_= nullptr;
+
+public:
+
+    TypeInfo("postprocessData");
+
+    /// @brief  Construct from systemControl and a boolean flag
+    ///  this constructor is used when postprocesing is active
+    ///  during simulation.
+    /// @param control const reference to systemControl 
+    postprocessData(const systemControl& control);
+
+    ~postprocessData()override = default;
+
+    add_vCtor
+    (
+        auxFunctions,
+        postprocessData,
+        systemControl
+    );
+
+    bool execute() override;
+    
+    
+     
+    bool wirte()const override
+    {
+        return false;
+    }
+
+};
+
+} // namespace pFlow
+
+#endif // __postprocessData_hpp__

src/PostProcessData/postprocessDataDict

@@ -0,0 +1,108 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName 	processDataDict;
+objectType 	dictionary;;
+fileFormat  ASCII;
+/*---------------------------------------------------------------------------*/
+
+runTimeActive 		yes;
+
+defaultTimeControl
+{
+	timeControl;
+	startTime;
+	endTime;
+	actionInterval  0.05; 
+}
+
+components
+(
+    velocityProb
+    {
+        method 	particleProbe;
+        region  idSelecttion;
+        field   velocity;
+        ids (0 10 100);
+        timeControl 	timeStep;
+        startTime 		0;
+        endTime 		infinity;
+        probInterval 	1;
+    }
+
+    comp2
+    {
+        method    uniformDistribution;
+        region spehre;
+        
+        sphereInfo
+        {
+            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
+    {
+        
+        region  line;
+        lineInfo
+        {
+            p1 ();
+            p2 ();
+            numPoints 10;
+            radius 0.01;
+        }
+        timeControl settingsDict; //default; 
+        type    numberBased;
+        operations();
+    }
+
+    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;
+                
+            }
+        );
+    }
+);

src/PostProcessData/postprocessTimeControl.hpp

@@ -0,0 +1,65 @@
+/*------------------------------- 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 __postprocessTimeControl_hpp__
+#define __postprocessTimeControl_hpp__
+
+#include "baseTimeControl.hpp"  
+
+namespace pFlow
+{
+class postprocessTimeControl 
+: 
+    public baseTimeControl
+{
+public:
+postprocessTimeControl(
+        const dictionary&          dict,
+        const baseTimeControl&     defaultTimeControl,
+        const baseTimeControl&     settingTimeControl
+    )
+    : 
+        baseTimeControl(0, 1, 1)
+    {
+        auto tControl = dict.getValOrSet<word>("timeControl", "default");
+        if(tControl == "default")
+        {
+            baseTimeControl::operator=(
+                defaultTimeControl
+            );
+        }
+        else if(tControl == "settingsDict")
+        {
+            baseTimeControl::operator=(
+                settingTimeControl
+            );
+        }
+        else
+        {
+            baseTimeControl::operator=( baseTimeControl(dict, "execution") );
+        }
+    }
+
+    // Additional methods and members can be added here
+};
+
+} // namespace pFlow
+
+#endif // __postprocessTimeControl_hpp__

src/PostProcessData/region/regionFields.hpp/regionField.hpp

@@ -0,0 +1,95 @@
+/*------------------------------- 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 __regionField_hpp__
+#define __regionField_hpp__
+
+#include "types.hpp"
+#include "regionPoints.hpp"
+#include "Field.hpp"
+
+namespace pFlow
+{
+
+template<typename T> 
+class regionField
+{
+private:
+           
+    /// the field value
+    Field<T, HostSpace>      field_;
+    
+    /// the region points
+    const regionPoints& regionPoints_;
+public:
+
+    TypeInfoTemplateNV11("regionField", T);
+
+    regionField(
+        const word& name, 
+        const regionPoints& rPoints,
+        const T defaultVal);
+    
+    regionField(const regionField&) = default;
+
+    regionField(regionField&&) = default;
+
+    regionField& operator=(const regionField&) = default;
+
+    regionField& operator=(regionField&&)  = default;
+
+    ~regionField() = default;
+
+    /// get the field value
+    T& operator[] (const uint32 i)
+    {
+        return field_[i];
+    }
+
+    /// get the field value
+    const T& operator[] (const uint32 i)const
+    {
+        return field_[i];
+    }
+
+    /// get field name 
+    word name()const
+    {
+        return field_.name();
+    }
+
+    /// @brief  write the field to a file 
+    /// @param os output file stream
+    /// @return true if successful and false if fails 
+    bool writeFieldToFile(iOstream& os)const;
+
+    /// @brief  write the field to vtk format (cell based)
+    /// @param os output file stream
+    /// @return true if successful and false if fails 
+    bool writeFieldToVtk(iOstream& os)const;
+
+};
+
+
+} // namespace pFlow
+
+#include "regionFieldTemplate.cpp"
+
+#endif // __regionField_hpp__

src/PostProcessData/region/regionFields.hpp/regionFieldTemplate.cpp

@@ -0,0 +1,10 @@
+
+template<typename T> 
+pFlow::regionField<T>::regionField(
+    const word& name, 
+    const regionPoints& rPoints,
+    const T defaultVal)
+:
+    field_(name, "regionFieldValue", rPoints.size(), rPoints.size(), defaultVal),
+    regionPoints_(rPoints)
+{}

src/PostProcessData/region/regionPoints/centerPointsRegionPoints.cpp

@@ -0,0 +1,86 @@
+#include "centerPointsRegionPoints.hpp"
+#include "fieldsDataBase.hpp"
+#include "Set.hpp"
+#include "pStructSelector.hpp"
+
+bool pFlow::centerPointsRegionPoints::selectIds()
+{
+    if(!firstTimeUpdate_) return true;
+    firstTimeUpdate_ = false;
+    
+    word selector =  probDict_.getVal<word>("selector");
+    
+    if(selector == "id")
+    {
+        auto idList = probDict_.getVal<uint32List>("ids");
+        Set<uint32> uniqueIds;
+
+        uniqueIds.insert(idList.begin(), idList.end());
+
+        for(auto& id:uniqueIds)
+        {
+            ids_.push_back(id);
+        }
+    }
+    else
+    {
+        auto selectorPtr = pStructSelector::create(
+            selector, 
+            database().pStruct(), 
+            probDict_);
+        auto selectedPoints = selectorPtr->selectedPoints();
+        ids_.resize(selectedPoints.size());
+        ids_.assign(selectedPoints.begin(), selectedPoints.end());
+    }
+
+    volume_.resize(ids_.size(),1.0);
+    diameter_.resize(ids_.size(),1.0);
+    center_.resize(ids_.size(), realx3(0,0,0));
+    selectedPoints_.resize(ids_.size(), -1);
+    return false;
+}
+
+pFlow::centerPointsRegionPoints::centerPointsRegionPoints(
+    const dictionary &dict,
+    fieldsDataBase &fieldsDataBase)
+    : regionPoints(dict, fieldsDataBase),
+      idName_(dict.getValOrSet<word>("idName", "id")),
+      probDict_(dict)
+{
+    auto idList = dict.getVal<uint32List>("ids");
+    Set<uint32> uniqueIds;
+
+    uniqueIds.insert(idList.begin(), idList.end());
+
+    for(auto& id:uniqueIds)
+    {
+        ids_.push_back(id);
+    }
+
+    selectedPoints_.resize(ids_.size());
+}
+
+bool pFlow::centerPointsRegionPoints::update()
+{
+    if(!selectIds()) return false;
+    
+    if(ids_.empty()) return true;
+
+    const auto& idField = database().updateFieldUint32(idName_);
+    selectedPoints_.fill(-1);
+
+    for(uint32 i = 0; i < idField.size(); ++i)
+    {
+        for( uint32 j=0; j< ids_.size(); ++j)
+        {
+            if(idField[i] == ids_[j])
+            {
+                selectedPoints_[j] = i;
+                break;
+            }
+        }
+    }
+
+    return false;
+}
+    

src/PostProcessData/region/regionPoints/centerPointsRegionPoints.hpp

@@ -0,0 +1,116 @@
+/*------------------------------- 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 __centerPointsRegionPoints_hpp__
+#define __centerPointsRegionPoints_hpp__
+
+#include "regionPoints.hpp"
+
+namespace pFlow
+{
+
+class centerPointsRegionPoints
+:
+    public regionPoints
+{
+private:
+    
+    bool            firstTimeUpdate_ = true;
+
+    /// the ids provided in the dictionary
+    uint32Vector    ids_;   
+
+    /// the volume of region
+    realVector       volume_ {"volume"};
+
+    realVector       diameter_{"diameter"};
+    /// center point of the region
+    /// this is not used in the idSelection region
+    realx3Vector     center_ {"center"};
+
+    /// the point indices that are selected by this region
+    uint32Vector     selectedPoints_{"selectedPoints"};
+    
+    /// the name of the id field
+    word             idName_= "id";
+
+    /// keeps the dictionary for first update use 
+    dictionary       probDict_;
+    
+    
+    bool selectIds();
+
+public:
+
+    TypeInfo("centerPoints");
+
+    centerPointsRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    ~centerPointsRegionPoints() override = default;
+
+    uint32 size()const override
+    {
+        return selectedPoints_.size();
+    }
+
+    bool empty()const override
+    {
+        return selectedPoints_.empty();
+    }
+
+    span<const real> volumes()const override
+    {
+        return span<const real>(volume_.data(), volume_.size());
+    }
+
+    span<const real> eqDiameters()const override
+    {
+        return span<const real>(diameter_.data(), diameter_.size());
+    }
+
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(center_.data(), center_.size());
+    }
+
+    span<const uint32> indices(uint32 elem)const override
+    {
+        return span<const uint32>(selectedPoints_.data(), selectedPoints_.size());
+    }
+
+    span<uint32> indices(uint32 elem) override
+    {
+        return span<uint32>(selectedPoints_.data(), selectedPoints_.size());
+    }
+
+    /// @brief  update the selected points based on the ids
+    /// @return true if the operation is successful
+    bool update() override;
+
+
+}; // class centerPointsRegionPoints
+
+} // namespace pFlow
+
+
+
+#endif // __centerPointsRegionPoints_hpp__

src/PostProcessData/region/regionPoints/lineRegionPoints.cpp

@@ -0,0 +1,83 @@
+#include "lineRegionPoints.hpp"
+#include "fieldsDataBase.hpp"
+
+pFlow::lineRegionPoints::lineRegionPoints
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDataBase
+)
+:
+    regionPoints(dict, fieldsDataBase),
+    line_(dict.subDict("lineInfo")),
+    centerPoints_("centerPoints"),
+    volumes_("volumes"),
+    selectedPoints_("selectedPoints")
+{
+    const auto& lDict = dict.subDict("lineInfo");
+    uint32 nPoints = lDict.getValMax<uint32>("numPoints",2);
+    realList raddi;
+
+    if( lDict.containsDataEntry("radii"))
+    {
+        raddi = lDict.getVal<realList>("radii");
+    }
+    else
+    {
+        auto r = lDict.getVal<real>("radius");
+        raddi = realList(nPoints, r);
+    }
+
+    if(raddi.size() != nPoints)
+    {
+        fatalErrorInFunction
+            << "The number elements in of radii list should be equal to the number of points"<<endl;
+            fatalExit;
+    }
+    
+    sphereRegions_.resize(nPoints, sphere(realx3(0,0,0),1));
+    centerPoints_.resize(nPoints);
+    volumes_.resize(nPoints);
+    selectedPoints_.resize(nPoints);
+    real dt = 1.0/(nPoints-1);
+    for(uint32 i = 0; i < nPoints; ++i)
+    {
+        centerPoints_[i] = line_.point(i*dt);
+        sphereRegions_[i] = pFlow::sphere(centerPoints_[i], raddi[i]);
+        volumes_[i] = sphereRegions_[i].volume();
+        diameters_[i] = 2*sphereRegions_[i].radius();
+    }
+}
+
+pFlow::span<const pFlow::uint32> pFlow::lineRegionPoints::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::lineRegionPoints::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;
+}

src/PostProcessData/region/regionPoints/lineRegionPoints.hpp

@@ -0,0 +1,100 @@
+/*------------------------------- 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 __lineRegionPoints_hpp__
+#define __lineRegionPoints_hpp__
+
+#include "regionPoints.hpp"
+#include "line.hpp"
+#include "sphere.hpp"
+#include "Vectors.hpp"
+
+namespace pFlow
+{
+
+class lineRegionPoints
+:
+    public regionPoints
+{
+private:
+
+    /// line region for selecting points
+    line                line_;
+
+    /// all sphere regions 
+    Vector<sphere>        sphereRegions_;
+
+    /// center poitns of regions/elements
+    realx3Vector          centerPoints_;
+
+    /// volumes of all elements/regions
+    realVector            volumes_;
+
+    realVector            diameters_;
+
+    /// the point indices that are selected by this region
+    Vector<uint32Vector>  selectedPoints_;
+
+public:
+
+    TypeInfo(line::TYPENAME());
+
+    lineRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    ~lineRegionPoints() override = default;
+
+    uint32 size()const override
+    {
+        return sphereRegions_.size();
+    }
+
+    bool empty()const override
+    {
+        return sphereRegions_.empty();
+    }
+
+    span<const real> volumes()const override
+    {
+        return span<const real>(volumes_.data(), volumes_.size());
+    }
+
+    span<const real> eqDiameters()const override
+    {
+        return span<const real>(diameters_.data(), diameters_.size());
+    }
+
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(centerPoints_.data(), centerPoints_.size());
+    }
+
+    span<const uint32> indices(uint32 elem)const override;
+    
+
+    bool update() override;
+
+};
+    
+}
+
+
+#endif // __lineRegionPoints_hpp__

src/PostProcessData/region/regionPoints/regionPoints.cpp

@@ -0,0 +1,28 @@
+#include "regionPoints.hpp"
+#include "fieldsDataBase.hpp"
+#include "Time.hpp"
+
+pFlow::regionPoints::regionPoints
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDataBase
+)
+:
+    fieldsDataBase_(fieldsDataBase)
+{}
+
+const pFlow::Time& pFlow::regionPoints::time() const
+{
+    return fieldsDataBase_.time();
+}
+
+const pFlow::fieldsDataBase & pFlow::regionPoints::database() const
+{
+    return fieldsDataBase_;
+}
+
+pFlow::fieldsDataBase& pFlow::regionPoints::database()
+{
+    return fieldsDataBase_;
+}
+

src/PostProcessData/region/regionPoints/regionPoints.hpp

@@ -0,0 +1,103 @@
+/*------------------------------- 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 __regionPoints_hpp__
+#define __regionPoints_hpp__
+
+#include "dictionary.hpp"
+#include "pointStructure.hpp"
+
+
+namespace pFlow
+{
+
+class fieldsDataBase;
+class Time;
+
+class regionPoints
+{
+    using PointsTypeHost = typename pointStructure::PointsTypeHost;
+
+    fieldsDataBase& fieldsDataBase_;
+
+public:
+
+    TypeInfo("regionPoints");
+
+
+    regionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+    
+    virtual ~regionPoints() = default;
+
+    const Time& time()const;
+
+    const fieldsDataBase& database()const;
+
+    fieldsDataBase& database();
+    
+    /// @brief size of elements 
+    virtual 
+    uint32 size()const = 0;
+    
+    /// @brief  check if the region is empty
+    virtual 
+    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 
+    virtual 
+    span<const real> volumes()const =0;
+    
+    virtual
+    span<const real> eqDiameters()const = 0; 
+
+    /// center points of elements 
+    virtual
+    span<const realx3> centers()const = 0;
+
+    /// indices of particles inside the element @var elem
+    virtual 
+    span<const uint32> indices(uint32 elem)const = 0;
+
+    virtual
+    span<uint32> indices(uint32 elem) = 0;
+
+    virtual 
+    bool update() = 0;
+    
+    /*virtual 
+    bool write()const=0;*/
+
+    /*static
+    uniquePtr<regionPoints> create(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);*/
+    
+};
+
+}
+
+#endif // __regionPoints_hpp__

src/PostProcessData/region/regionPoints/sphereRegionPoints.cpp

@@ -0,0 +1,31 @@
+#include "sphereRegionPoints.hpp"
+#include "fieldsDataBase.hpp"
+
+pFlow::sphereRegionPoints::sphereRegionPoints
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDataBase
+)
+:
+    regionPoints(dict, fieldsDataBase),
+    sphereRegion_(dict.subDict("sphereInfo")),
+    volume_(sphereRegion_.volume()),
+    diameter_(2*sphereRegion_.radius()),
+    selectedPoints_("selectedPoints")
+{
+}
+
+bool pFlow::sphereRegionPoints::update()
+{
+    const auto points = database().updatePoints();
+    selectedPoints_.clear();
+    for(uint32 i = 0; i < points.size(); ++i)
+    {
+        if( sphereRegion_.isInside(points[i]))
+        {
+            selectedPoints_.push_back(i);
+        }
+    }
+ 
+    return false;
+}

src/PostProcessData/region/regionPoints/sphereRegionPoints.hpp

@@ -0,0 +1,98 @@
+/*------------------------------- 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 __sphereRegionPoints_hpp__
+#define __sphereRegionPoints_hpp__
+
+#include "regionPoints.hpp"
+#include "sphere.hpp"
+#include "Vectors.hpp"
+
+namespace pFlow
+{
+class sphereRegionPoints
+:
+    public regionPoints
+{
+private:
+
+    /// spehre region for selecting points
+    sphere  sphereRegion_;
+
+    /// the volume of region
+    real    volume_;
+
+    real    diameter_;
+
+    /// the point indices that are selected by this region
+    uint32Vector    selectedPoints_;
+
+public:
+
+    TypeInfo(sphere::TYPENAME());
+
+    sphereRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    ~sphereRegionPoints() override = default;
+    
+    uint32 size()const override
+    {
+        return 1;
+    }
+
+    bool empty()const override
+    {
+        return false;
+    }
+
+    span<const real> volumes()const override
+    {
+        return span<const real>(&volume_, 1);
+    }
+
+    span<const real> eqDiameters()const override
+    {
+        return span<const real>(&diameter_, 1);
+    }
+
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(&sphereRegion_.center(), 1);
+    }
+        
+    span<const uint32> indices(uint32 elem)const override
+    {
+        return span<const uint32>(selectedPoints_.data(), selectedPoints_.size());
+    }
+
+    span<uint32> indices(uint32 elem) override
+    {
+        return span<uint32>(selectedPoints_.data(), selectedPoints_.size());
+    }
+
+    bool update()override;
+    
+};
+
+}
+
+#endif // __sphereRegionPoints_hpp__