Merge branch 'PhasicFlow:main' into benchmarks

cmake/autoComplete

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

src/CMakeLists.txt

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

src/Integration/AdamsBashforth2/AdamsBashforth2.cpp

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

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

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

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

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

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

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

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

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

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

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

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

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

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

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

src/Integration/integration/integration.cpp

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

src/Integration/integration/integration.hpp

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

src/Particles/GrainParticles/grainParticles/grainParticles.cpp

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

src/Particles/GrainParticles/grainShape/grainShape.cpp

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

src/Particles/GrainParticles/grainShape/grainShape.hpp

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

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

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

src/Particles/SphereParticles/sphereShape/sphereShape.cpp

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

src/Particles/SphereParticles/sphereShape/sphereShape.hpp

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

src/Particles/dynamicPointStructure/dynamicPointStructure.cpp

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

src/Particles/particles/particles.hpp

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

src/Particles/particles/shape/shape.cpp

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

src/Particles/particles/shape/shape.hpp

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

src/PostprocessData/CMakeLists.txt

@@ -0,0 +1,35 @@
+set(SourceFiles 
+    # Main postprocess data
+    postprocessData/postprocessData.cpp
+
+    # Fields database
+    fieldsDataBase/fieldsDataBase.cpp
+    fieldsDataBase/simulationFieldsDataBase.cpp
+  
+    # Regions
+    region/regionPoints/regionPoints/regionPoints.cpp
+    region/regionPoints/sphereRegionPoints/sphereRegionPoints.cpp
+    region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
+    region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
+    region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
+
+    # Postprocess components
+    postprocessComponent/postprocessComponent/postprocessComponent.cpp
+    postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp 
+    postprocessComponent/PostprocessComponent/PostprocessComponents.cpp
+
+    # Operations
+    operation/postprocessOperation/postprocessOperation.cpp
+    operation/PostprocessOperation/PostprocessOperationSum.cpp
+    operation/PostprocessOperation/PostprocessOperationAverage.cpp
+    operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp
+    
+    operation/includeMask/includeMask.cpp
+    operation/includeMask/IncludeMasks.cpp
+
+    
+)
+
+set(link_libs Kokkos::kokkos phasicFlow Particles)
+
+pFlow_add_library_install(PostprocessData SourceFiles link_libs)

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,855 @@
+/*------------------------------- 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 "systemControl.hpp"
+#include "fieldsDataBase.hpp"
+#include "fieldFunctions.hpp"
+
+
+namespace pFlow
+{
+
+bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+
+}
+
+bool pFlow::fieldsDataBase::loadPointStructureToTime()
+{
+    return false;
+}
+
+pFlow::word pFlow::fieldsDataBase::getPointFieldType(const word &name) const
+{
+    word pfType =  time_.lookupObjectTypeName(name);
+    word type, space;
+    if(!pointFieldGetType(pfType, type, space))
+    {
+        fatalErrorInFunction
+            <<"Error in retriving the type of pointField "
+            << pfType<<endl;
+        fatalExit;
+    }
+    return type;
+}
+
+bool pFlow::fieldsDataBase::checkForUpdate(const word &compoundName, bool forceUpdate)
+{
+    auto t = currentTime();
+    bool shouldUpdate = false;
+    
+    if(auto [iter, found]= captureTime_.findIf(compoundName); found)
+    {
+        shouldUpdate = iter->second < t || forceUpdate;
+        iter->second = t;
+    }
+    else
+    {
+        shouldUpdate = true;
+        captureTime_.insertIf(compoundName, 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());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetVolume(bool forceUpdate)
+{
+    const word fName = "volume";
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+        
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        auto vols = getShape().volume();
+        
+        FieldTypeHost<real> volField
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< volField.size(); i++)
+        {
+            volField[i] = vols[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(volField)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetDensity(bool forceUpdate)
+{
+    const word fName = "density";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        const auto dens = getShape().density();
+
+        FieldTypeHost<real> denFeild
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< denFeild.size(); i++)
+        {
+            denFeild[i] = dens[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(denFeild)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetOne(bool forceUpdate)
+{
+    const word fName = "one";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            FieldTypeHost<real> 
+            (
+                fName, 
+                "value",
+                pointFieldSize(),
+                1.0
+            )  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+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(systemControl& control, bool inSimulation)
+:
+    time_(control.time()),
+    inSimulation_(inSimulation)
+{}
+
+pFlow::timeValue pFlow::fieldsDataBase::currentTime() const
+{
+    return time_.currentTime();
+}
+
+bool pFlow::fieldsDataBase::getFieldTypeNameFunction
+(
+    const word& compoundName,
+    word& pointFieldName, 
+    word& originalType, 
+    word& typeAfterFunction,
+    Functions& func
+)const
+{
+    if( !findFunction(compoundName, pointFieldName, func))
+    {
+        fatalErrorInFunction
+            <<"Error in processing function for field: "
+            <<compoundName<<endl;;
+        return false;
+    }
+
+    if( reservedFieldNames_.contains(pointFieldName))
+    {
+        // The name is in the reserved fields 
+        originalType = reservedFieldNames_.find(pointFieldName)->second;
+    }
+    else
+    {
+        // the name is in the Time object 
+        if(pointFieldNameExists(pointFieldName))
+        {
+            originalType = getPointFieldType(pointFieldName);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "The field name: "<< pointFieldName
+                << " is not in the Time object.\n"
+                <<" Avaiable ones are: \n"<< time().objectNames()<<endl;
+            return false;
+        }
+    }
+
+    if(!inputOutputType(func, originalType, typeAfterFunction))
+    {
+        fatalErrorInFunction
+            <<"Cannnot determine the input and output types for: "
+            << compoundName<<endl;
+        return false;
+    }
+
+    return true;
+}
+
+bool pFlow::fieldsDataBase::getFieldType
+(
+    const word& compoundName,
+    word& originalType, 
+    word& typeAfterFunction
+) const
+{
+    Functions func;
+    word fieldName;
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
+}
+
+bool pFlow::fieldsDataBase::getFieldType
+(
+    const word &compoundName, 
+    word &typeAfterFunction
+) const
+{
+    Functions func;
+    word originalType;
+    word fieldName;
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
+}
+
+pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updatePoints(bool forceUpdate)
+{
+    const word fName = "position";
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+        
+    if(shouldUpdate)
+    {
+        // load pointStructure 
+        if(!loadPointStructureToTime())
+        {
+            fatalErrorInFunction<< "Error in loading pointStructure to Time object."<<endl;
+            fatalExit;
+        }
+        const auto& pstruct = pStruct(); 
+        allFields_.emplaceBackOrReplace<PointsTypeHost>
+        (
+            fName, 
+            pstruct.activePointsHost()
+        );
+    }
+   
+    auto& points = allFields_.getObject<PointsTypeHost>(fName);
+    
+    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( !getFieldTypeNameFunction(
+        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( !getFieldTypeNameFunction(
+        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( !getFieldTypeNameFunction(
+        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>())
+    {
+        if(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 typeAfterFunction, originalType;
+    
+    if( !getFieldType(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 updateFieldRealx3(compoundName, forceUpdate);
+    }
+    else if( typeAfterFunction == getTypeName<realx4>() )
+    {
+        return updateFieldRealx4(compoundName, forceUpdate);
+    }
+    else if( typeAfterFunction == getTypeName<real>() )
+    {
+        return updateFieldReal(compoundName, forceUpdate);
+    }
+    else
+    {
+        fatalErrorInFunction<< "Invalid feild "<< compoundName<<endl;
+        fatalExit;
+        return span<real>(nullptr, 0);
+    }
+}
+
+
+
+pFlow::uniquePtr<pFlow::fieldsDataBase> 
+    pFlow::fieldsDataBase::create(systemControl& control, bool inSimulation)
+{
+    word dbType;
+    if(inSimulation)
+    {
+        dbType = "fieldsDataBase<simulation>";
+    }
+    else
+    {
+        dbType = "fieldsDataBase<postSimulation>";
+    }
+    
+    if( boolvCtorSelector_.search(dbType) )
+    {
+        auto objPtr = 
+            boolvCtorSelector_[dbType](control, inSimulation);
+        return objPtr;
+    }
+    else
+    {
+        printKeys
+        ( 
+            fatalError << "Ctor Selector "<< 
+            dbType << " does not exist. \n"
+            <<"Available ones are: \n\n"
+            ,
+            boolvCtorSelector_
+        );
+        fatalExit;
+    }
+    
+    return nullptr;
+}
+
+bool pFlow::pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace)
+{
+	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,277 @@
+/*------------------------------- 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 "fieldsDataBaseDclr.hpp"
+
+#include "pointStructure.hpp"
+#include "pointFields.hpp"
+#include "anyList.hpp"
+#include "Map.hpp"
+#include "shape.hpp"
+
+
+namespace pFlow
+{
+
+class systemControl;
+class Time;
+
+
+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
+
+        /// List to store all the point fields 
+        anyList             allFields_;
+
+        /// Map to store the last capture time of each field
+        wordMap<timeValue>  captureTime_;
+
+        /// Reference to the Time object
+        Time&               time_; 
+
+        /// Flag indicating if we're in simulation mode
+        bool                inSimulation_ = false;
+
+protected:
+
+        /// Map of reserved field names with their corresponding types
+        static const inline std::map<word, word> reservedFieldNames_
+        {
+            {"position", "realx3"},
+            {"one", "real"},
+            {"volume", "real"},
+            {"density", "real"}
+        };
+
+        /// check if pointField name exists in Time or time folder 
+        virtual 
+        bool pointFieldNameExists(const word& name)const = 0;
+        
+        /// Loads a pointField with a given name to the Time object 
+        virtual 
+        bool loadPointFieldToTime(const word& name)= 0;
+
+        virtual 
+        bool loadPointStructureToTime()=0;
+
+        /// get the type name of the pointField in the Time object 
+        word getPointFieldType(const word& name)const;
+
+        /// Checks if a field needs to be updated based on capture time
+        bool checkForUpdate(const word& compoundName, bool forceUpdate = false);
+
+        template<ValidFieldType T>
+        span<T> updateField(const word& name, bool forceUpdate = false);
+
+        template<ValidFieldType T>
+        span<T> updateReservedField(const word& name, bool forceUpdate = false);
+        
+        span<real> createOrGetRealField(const word& name);
+
+        span<real> createOrGetVolume(bool forceUpdate=false);
+
+        span<real> createOrGetDensity(bool forceUpdate=false);
+
+        span<real> createOrGetOne(bool forceUpdate=false);
+
+        static 
+        bool findFunction(
+            const word& compoundFieldName, 
+            word& fieldName, 
+            fieldsDataBase::Functions& func );
+        
+        static
+        bool inputOutputType(
+            fieldsDataBase::Functions func, 
+            const word& inputType, 
+            word& outputType);
+
+
+protected:
+    
+    // - protected member functions
+    
+
+    virtual 
+    bool checkTimeFolder(const word& fieldName) const = 0;
+    
+    virtual 
+    const shape& getShape() const= 0;
+
+
+    /// @brief return the size of pointStructure 
+    uint32 pointFieldSize()
+    {
+        auto s = updatePoints();
+        return s.size();
+    }
+
+public:
+    
+    // - Type info
+            
+        TypeInfo("fieldsDataBase");
+
+    // - constructors 
+
+        fieldsDataBase(systemControl& control, 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
+        virtual ~fieldsDataBase() = default;
+
+        create_vCtor
+        (
+            fieldsDataBase,
+            bool,
+            (systemControl& control, bool inSimulation),
+            (control, inSimulation)
+        );
+  
+    
+    // - Public Access Functions
+        /// returns the current time
+        timeValue currentTime()const;
+
+        /// const ref to object Time 
+        const Time& time()const 
+        {
+            return time_;
+        }
+        
+        Time& time()
+        {
+            return time_;
+        }
+
+    // - Public Member Functions
+
+        bool getFieldTypeNameFunction
+        (
+            const word& compoundName,
+            word& pointFieldName, 
+            word& originalType, 
+            word& typeAfterFunction,
+            Functions& func
+        )const;
+
+        
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& originalType, 
+            word& typeAfterFunction
+        ) const;
+
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& typeAfterFunction
+        ) const;
+
+        
+        /// update pointStructure if necessary 
+        span<realx3> updatePoints(bool forceUpdate = false);
+        
+        /// update a field with realx3 type 
+        span<realx3> updateFieldRealx3(
+            const word& compoundName, 
+            bool forceUpdate = false);
+        
+        span<realx4> updateFieldRealx4(
+            const word& compoundName, 
+            bool forceUpdate = false);
+        
+        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);
+        
+        virtual 
+        const pointStructure& pStruct()const = 0;
+
+        virtual 
+        void resetTimeFolder() = 0;
+    
+    static
+    uniquePtr<fieldsDataBase> create(
+        systemControl& control, 
+        bool inSimulation);
+};
+
+} // namespace pFlow
+
+#include "fieldsDataBaseTemplates.cpp"
+
+#endif //__fieldsDataBased_hpp__

src/PostprocessData/fieldsDataBase/fieldsDataBaseDclr.hpp

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

src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp

@@ -0,0 +1,141 @@
+/*------------------------------- 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(reservedFieldNames_.contains(name))
+        {
+            return updateReservedField<T>(name, true);
+        }   
+        else
+        {
+            if( loadPointFieldToTime(name) )
+            {
+                const auto& pField = time_.lookupObject<pointField_D<T>>(name);
+                allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
+                name, 
+                pField.activeValuesHost());
+            }
+            else
+            {
+                fatalErrorInFunction<<"Error in loading the pointField "<<name<<endl;
+                fatalExit;
+            }
+            
+        }
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<T>>(name);
+    
+    return span<T>(
+        field.data(), 
+        field.size());
+    
+}
+
+
+template<pFlow::ValidFieldType T>
+inline
+pFlow::span<T> pFlow::fieldsDataBase::updateReservedField
+(
+    const word& name, 
+    bool forceUpdate
+)
+{
+    if(name == "one")
+    {   
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetOne(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field one."<<endl;
+            fatalExit;
+        }
+    }
+    else if(name == "volume")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetVolume(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field volume."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "density")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetDensity(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field density."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "position")
+    {
+        if constexpr( std::same_as<T, realx3>)
+        {
+            return updatePoints(forceUpdate);            
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field position."<<endl;
+            fatalExit;
+        }
+    }
+    
+    fatalErrorInFunction<<"Not supported field name "<<endl;
+    fatalExit;
+    return span<T>(nullptr, 0);
+
+}
+
+#endif //__fieldsDataBaseTemplates_hpp__

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.cpp

@@ -0,0 +1,62 @@
+#include "Time.hpp"
+#include "simulationFieldsDataBase.hpp"
+#include "dynamicPointStructure.hpp"
+#include "vocabs.hpp"
+
+namespace pFlow
+{
+    bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+}
+
+bool pFlow::simulationFieldsDataBase::pointFieldNameExists(const word &name) const
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::simulationFieldsDataBase::loadPointFieldToTime(const word &name)
+{
+    return time().lookupObjectName(name);
+}
+
+bool pFlow::simulationFieldsDataBase::loadPointStructureToTime()
+{
+    // it is already in the Time object 
+    return time().lookupObjectName(pointStructureFile__);
+}
+
+bool pFlow::simulationFieldsDataBase::checkTimeFolder(const word &fieldName) const
+{
+    return true;
+}
+
+const pFlow::shape& pFlow::simulationFieldsDataBase::getShape() const
+{
+    return shape_;
+}
+
+pFlow::simulationFieldsDataBase::simulationFieldsDataBase
+(
+    systemControl &control, 
+    bool inSimulation
+)
+:
+    fieldsDataBase(control, inSimulation),
+    shape_
+    (
+        dynamic_cast<const shape&>(*control.caseSetup().lookupObjectPtr(shapeFile__))
+    )
+{
+}
+
+const pFlow::pointStructure &pFlow::simulationFieldsDataBase::pStruct() const
+{
+    return 
+    static_cast<const pointStructure&>
+    (
+        time().lookupObject<dynamicPointStructure>(pointStructureFile__)
+    );
+}
+
+void pFlow::simulationFieldsDataBase::resetTimeFolder()
+{
+}

src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.hpp

@@ -0,0 +1,79 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __simulationFieldsDataBase_hpp__
+#define __simulationFieldsDataBase_hpp__
+
+#include "fieldsDataBase.hpp"
+
+namespace pFlow
+{
+
+class simulationFieldsDataBase 
+: 
+    public fieldsDataBase
+{
+private:
+
+    const   shape&       shape_;
+
+protected:
+    
+    /// check if pointField name exists in Time or time folder 
+    bool pointFieldNameExists(const word& name)const override;
+
+    /// Loads a pointField with a given name to the Time object. 
+    /// For simulation, it just checks if the name exists 
+    bool loadPointFieldToTime(const word& name) override;
+
+    /// Loads pointStructure to the Time object
+    /// For simulation, it just checks if the name exists
+    bool loadPointStructureToTime() override;
+
+    bool checkTimeFolder(const word& fieldName) const override;
+
+    const shape& getShape() const override;  
+    
+public:
+
+    TypeInfo("fieldsDataBase<simulation>");
+
+    simulationFieldsDataBase(systemControl& control, bool inSimulation);
+
+    ~simulationFieldsDataBase() override = default;
+
+    add_vCtor
+    (
+        fieldsDataBase,
+        simulationFieldsDataBase,
+        bool
+    );
+
+    const pointStructure& pStruct()const override;
+
+    void resetTimeFolder() override;
+
+};
+
+
+}
+
+
+#endif //__simulationFieldsDataBase_hpp__

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp

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

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp

@@ -0,0 +1,74 @@
+#include "PostprocessOperationSum.hpp"
+#include "dictionary.hpp"
+#include "fieldsDataBase.hpp"
+#include "operationFunctions.hpp"
+
+/// Constructs sum processor and initializes result field based on input field type
+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;
+    }
+}
+
+/// Performs weighted sum of field values within each region
+bool pFlow::PostprocessOperationSum::execute
+(
+    const std::vector<span<real>>& weights,
+    const regionField<real>& volFactor
+)
+{
+    auto allField = database().updateFieldAll(fieldName());
+    auto phi = database().updateFieldReal(
+        phiFieldName());
+
+    auto mask = getMask();
+    word procName = processedFieldName();
+    const auto& regP = regPoints();
+    bool dbVol = divideByVolume();
+
+    processedRegField_ = makeUnique<processedRegFieldType>
+    (   
+        std::visit([&](auto&& field)->processedRegFieldType
+            { 
+                return executeSumOperation(
+                procName,
+                field,
+                volFactor,
+                dbVol,
+                weights,
+                phi,
+                mask);
+            },
+            allField)
+    );
+   
+    return true;
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp

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

src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp

@@ -0,0 +1,195 @@
+/*------------------------------- 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 __operationFunctions_hpp__
+#define __operationFunctions_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 regionField<real>&        volFactor,
+    const bool                      devideByVol,
+    const std::vector<span<real>>&  weights,
+    const span<real>&               phi,
+    const includeMask::Mask&        mask
+)
+{
+    const auto& regPoints = volFactor.regPoints(); 
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        
+        auto w = weights[reg];
+        T sum = T{};
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( index!= -1 && mask( index ))
+            {
+                sum += w[n] * field[index]* phi[index];
+            }
+            n++;
+        }
+        if(devideByVol)
+        {
+            processedField[reg] = sum/(volFactor[reg] * vols[reg]);
+        }
+        else
+        {
+            processedField[reg] = sum;
+        }
+        
+    }
+
+    return processedField;
+}
+
+template<typename T>
+regionField<T> executeAverageOperation
+(
+    const word&                    regFieldName,
+    const span<T>&                 field, 
+    const regionField<real>&       volFactor,
+    const bool                     devideByVol,
+    const std::vector<span<real>>& weights,
+    const span<real>&              phi,
+    const includeMask::Mask&       mask
+)
+{
+    
+    const auto& regPoints = volFactor.regPoints(); 
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto w = weights[reg];
+        T sumNum = T{};
+        real sumDen = 0;
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( index!= -1)
+            {
+                if( mask(index))
+                {
+                    sumNum += w[n] * field[index]* phi[index];
+                }
+                sumDen += w[n] * phi[index];
+            }
+            
+            n++;
+        }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / (volFactor[reg] * vols[reg]);
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
+        
+    }
+
+    return processedField;
+}
+
+
+template<typename T>
+regionField<T> executeFluctuation2Operation
+(
+    const word&                    regFieldName,
+    const span<T>&                 field,
+    const regionField<T>&          fieldAvg,
+    const regionField<real>&       volFactor, 
+    const bool                     devideByVol, 
+    const std::vector<span<real>>& weights,
+    const includeMask::Mask&       mask
+)
+{
+    const auto& regPoints = fieldAvg.regPoints();
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto vols = regPoints.volumes();
+
+    for(uint32 reg =0; reg<regPoints.size(); reg++)
+    {
+        auto partIndices = regPoints.indices(reg);
+        auto w = weights[reg];
+        auto vol = volFactor[reg] * vols[reg];
+        T avField{};
+        if(devideByVol)
+        {
+            avField = vol * fieldAvg[reg];
+        }
+        else
+        {
+            avField = fieldAvg[reg];
+        }
+        
+        T sumNum = T{};
+        real sumDen = 0;
+        uint n = 0;
+        for(auto index:partIndices)
+        {
+            if( index!= -1)
+            {
+                if( mask(index))
+                {
+                    sumNum += w[n] * pow( avField- field[index],static_cast<real>(2));
+                }
+                sumDen += w[n];
+            }
+            n++;
+        }
+        
+        if(devideByVol)
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue) / vol;
+        }
+        else
+        {
+            processedField[reg] = sumNum / max(sumDen, smallValue);
+        }
+        
+    }
+
+    return processedField;
+}
+
+} // namespace pFlow
+
+#endif //__operationFunctions_hpp__

src/PostprocessData/operation/includeMask/IncludeMask.hpp

@@ -0,0 +1,277 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+
+/**
+ * @brief A template class implementing includeMask for filtering data based on field values
+ * 
+ * The IncludeMask class creates boolean masks that identify which elements from a field
+ * satisfy a given condition. It applies an operator to each element of a field and
+ * generates a mask (vector of booleans) where true means the element satisfies the
+ * condition and should be included.
+ * 
+ * @tparam T The data type of the field (real, realx3, realx4)
+ * @tparam Operator The operation class that defines the condition (e.g., greaterThanOp)
+ * 
+ * The class has a specialized implementation for allOp operator which includes all elements.
+ * 
+ * Usage example in postprocessDataDict:
+ * ```
+ * // Create a dictionary with the required configuration for filtering
+ * 
+ * {
+ *     includeMask lessThan;
+ *
+ *     // Diameter of par1 is 0.003, so these settings
+ *     // will select only particles of type par1
+ *     lessThanInfo
+ *     {
+ *         field   diameter;
+ *         value   0.0031;
+ *     }
+ *    
+ * }
+ * 
+ * ```
+ */
+
+#ifndef __IncludeMask_hpp__
+#define __IncludeMask_hpp__
+
+#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:
+
+    /// Type alias for the mask container returned by getMask()
+    using Mask = typename includeMask::Mask;
+
+private:
+	
+    /// Boolean vector the filtering status of each element (true = include)
+    std::vector<bool>       mask_;
+
+    /// Comparison operator instance that evaluates the filtering condition
+    Operator                operator_;
+
+    /// Name of the field to apply filtering on
+    word                    fieldName_;
+    
+    /// Timestamp when mask was last updated (-1 indicates never updated)
+    timeValue               lastUpdated_ = -1;
+
+    /// Updates the mask based on current field values if needed, returns true if successful
+    bool updateMask()
+    {
+        timeValue t = database().currentTime();
+        
+        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;
+    }
+
+    /// Returns the name of the operator as a string (from operator's TYPENAME)
+    static
+    word operatorName()
+    {
+        return Operator::TYPENAME();
+    }
+
+public:
+
+	TypeInfoTemplate12("IncludeMask", T, Operator);
+
+    /// Constructs an IncludeMask using settings from dictionary and field database
+	IncludeMask(
+		const dictionary& dict, 
+        fieldsDataBase& feildsDB)
+	:
+		includeMask(dict, feildsDB),
+		operator_(dict.subDict(operatorName()+"Info")),
+        fieldName_(
+            dict.subDict
+            (
+                operatorName()+"Info"
+            ).getVal<word>("field"))
+	{}
+    
+    IncludeMask(
+        const word& type,
+		const dictionary& dict, 
+        fieldsDataBase& feildsDB)
+	:
+		includeMask(type, dict, feildsDB),
+		operator_(dict.subDict(operatorName()+"Info")),
+        fieldName_(
+            dict.subDict
+            (
+                operatorName()+"Info"
+            ).getVal<word>("field"))
+	{}
+
+    /// Add virtual constructor pattern for creating instances
+	add_vCtor
+    (
+		includeMask,
+		IncludeMask,
+		dictionary
+    );
+
+    add_vCtor
+    (
+		includeMask,
+		IncludeMask,
+		word
+    );
+
+    /// Returns the mask for filtering elements (updates the mask if necessary)
+    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);
+    }
+
+    IncludeMask(
+        const word& type,
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB)
+    :
+        includeMask(type, opDict, feildsDB)
+    {
+        span<realx3>  s = database().updatePoints();        
+        mask_.resize(s.size(), true);
+    }
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        dictionary
+    );
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        word
+    );
+
+    Mask getMask()override
+    {
+        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,152 @@
+/*------------------------------- 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::includeMask::includeMask
+(
+	const word &type, 
+	const dictionary &opDict, 
+	fieldsDataBase &fieldsDB
+)
+:
+	database_(fieldsDB)
+{
+}
+
+pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
+(
+    const dictionary& opDict, 
+    fieldsDataBase& fieldsDB
+)
+{
+	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( !fieldsDB.getFieldType(maskField, fieldType) )
+		{
+			fatalErrorInFunction<<"Error in retriving the type of field"
+				<< maskField <<" from dictionary "
+				<< maskDict.globalName()
+				<< endl;
+			fatalExit;
+			return nullptr;
+		}
+	}
+	else
+	{
+		fieldType = getTypeName<real>();
+	}
+	
+	word method = angleBracketsNames2("IncludeMask", fieldType, mask);
+	
+	if( dictionaryvCtorSelector_.search(method) )
+	{
+		auto objPtr = 
+			dictionaryvCtorSelector_[method]
+			(opDict, fieldsDB);
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< 
+			method << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+		return nullptr;
+	}
+	return nullptr;
+}
+
+pFlow::uniquePtr<pFlow::includeMask> 
+	pFlow::includeMask::create
+(
+	const word &type, 
+	const dictionary &opDict, 
+	fieldsDataBase &fieldsDB
+)
+{
+	word fieldType;
+	if( type != "all")
+	{
+		auto& maskDict = opDict.subDict(type+"Info");
+		word maskField = maskDict.getVal<word>("field");
+		
+		if( !fieldsDB.getFieldType(maskField, fieldType) )
+		{
+			fatalErrorInFunction<<"Error in retriving the type of field"
+				<< maskField <<" from dictionary "
+				<< maskDict.globalName()
+				<< endl;
+			fatalExit;
+			return nullptr;
+		}
+	}
+	else
+	{
+		fieldType = getTypeName<real>();
+	}
+	
+	word method = angleBracketsNames2("IncludeMask", fieldType, type);
+	
+	if( wordvCtorSelector_.search(method) )
+	{
+		auto objPtr = 
+			wordvCtorSelector_[method]
+			(type, opDict, fieldsDB);
+		return objPtr;
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< 
+			method << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		fatalExit;
+		return nullptr;
+	}
+	return nullptr;
+}

src/PostprocessData/operation/includeMask/includeMask.hpp

@@ -0,0 +1,195 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+/**
+ * @class includeMask
+* @brief Base class for creating inclusion masks for data filtering
+* 
+ * includeMask is an abstract base class for creating masks that filter data
+ * from a fieldsDataBase. Derived classes implement specific filtering criteria
+ * through the getMask() method which returns a Mask object - a wrapper around
+ * a vector of booleans that indicates which elements to include/exclude.
+ * 
+ * This class follows a factory pattern with create() methods that instantiate
+ * the appropriate derived mask type based on dictionary settings.
+ * 
+ * Derived classes can implement various filtering strategies such as:
+ * - Filtering by field values
+ * - Filtering by spatial regions
+ * - Filtering by predefined criteria
+ * 
+ * The Mask objects returned by getMask() can be used in postprocessing operations
+ * to include only the data points that match the specified criteria.
+ */
+
+#ifndef __includeMask_hpp__
+#define __includeMask_hpp__
+
+#include <vector> 
+
+#include "virtualConstructor.hpp"
+
+namespace pFlow
+{
+
+// forward declaration
+class fieldsDataBase;
+class dictionary;
+
+
+class includeMask
+{
+public:
+
+	/// @brief Wrapper around a boolean vector that represents elements to include/exclude
+	/// Provides a functional interface to access the underlying mask vector
+	/// where true indicates inclusion and false indicates exclusion.
+	class Mask
+	{
+		/// @brief Reference to the underlying boolean vector
+		const std::vector<bool>& mask_;
+
+		public:
+
+		/// @brief Constructor from a boolean vector
+		/// @param msk Boolean vector where true means include, false means exclude
+		Mask(const std::vector<bool>& msk)
+		:
+			mask_(msk)
+		{}
+
+		/// @brief Copy constructor
+		Mask(const Mask&) = default;
+
+		/// @brief Move constructor
+		Mask(Mask&&) = default;
+				
+		/// @brief Destructor
+		~Mask()=default;
+
+		/// @brief Get the size of the mask
+		auto size()const
+		{
+			return mask_.size();
+		}
+
+		/// @brief Check if element at index i should be included
+		/// @param i Index to check
+		bool operator()(uint32 i)const
+		{
+			return mask_[i];
+		}
+	};
+
+private:
+	
+	/// @brief Reference to the database containing fields to be filtered
+    fieldsDataBase&     database_;
+	
+public:
+
+	TypeInfo("includeMask");
+
+	/// @brief Constructor
+	/// @param opDict Dictionary containing operation settings
+	/// @param feildsDB Database of fields 
+	includeMask(const dictionary& opDict, fieldsDataBase& feildsDB);
+
+	/// @brief Constructor with explicit type
+	/// @param type Type of mask to create
+	/// @param opDict Dictionary containing operation settings
+	/// @param feildsDB Database of fields 
+    includeMask(const word& type, const dictionary& opDict, fieldsDataBase& feildsDB);
+
+	/// @brief Virtual destructor
+	virtual ~includeMask() = default;
+
+	/// @brief Virtual constructor pattern implementation using dictionary
+    create_vCtor
+    (
+        includeMask,
+        dictionary,
+        (
+            const dictionary& opDict, fieldsDataBase& feildsDB
+        ),
+        (opDict, feildsDB)		
+    );
+
+	/// @brief Virtual constructor pattern implementation using type and dictionary
+    create_vCtor
+    (
+        includeMask, 
+        word,
+        (
+            const word& type, 
+            const dictionary& opDict, 
+            fieldsDataBase& feildsDB
+        ),
+        (type, opDict, feildsDB)
+    );
+	 
+	/// @brief Get const access to the database
+	/// @return Const reference to the fields database
+	const fieldsDataBase& database()const
+	{
+		return database_;
+	}
+
+	/// @brief Get non-const access to the database
+	/// @return Reference to the fields database
+    fieldsDataBase& database()
+    {
+        return database_;
+    }
+
+	/// @brief Get the mask for filtering elements
+	/// @return Mask object indicating which elements to include
+    virtual 
+    Mask getMask()= 0;
+
+	/// @brief Factory method to create appropriate mask type from dictionary
+	/// @param opDict Dictionary with mask settings including type
+	/// @param feildsDB Database of fields to filter
+	/// @return Unique pointer to created mask instance
+	static 
+	uniquePtr<includeMask> create(
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB);
+    
+	/// @brief Factory method to create mask with explicit type
+	/// @param type Type of mask to create
+	/// @param opDict Dictionary with mask settings
+	/// @param feildsDB Database of fields to filter
+	/// @return Unique pointer to created mask instance
+    static
+    uniquePtr<includeMask> create(
+        const word& type,
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB);
+	
+};
+
+
+
+} // pFlow
+
+#endif //__IncludeMask_hpp__
+
+

src/PostprocessData/operation/includeMask/maskOperations.hpp

@@ -0,0 +1,181 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __maskOperation_hpp__
+#define __maskOperation_hpp__
+
+#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/postprocessOperation.cpp

@@ -0,0 +1,211 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "Time.hpp"
+#include "postprocessOperation.hpp"
+#include "regionPoints.hpp"
+#include "fieldsDataBase.hpp"
+
+namespace pFlow
+{
+
+template<typename T>
+inline
+bool writeField
+(
+    iOstream&             os,
+    timeValue             t,
+    const regionField<T>  field,
+    uint32                threshold,
+    const T&              defValue=T{}    
+)
+{
+    const auto& regPoints = field.regPoints();
+    const uint32 n = field.size();
+    os<<t<<tab;
+    for(uint32 i=0; i<n; i++)
+    {
+        auto numPar = regPoints.indices(i).size();
+        if(numPar >= threshold)
+        {
+            if constexpr(std::is_same_v<T,realx3>)
+            {
+                os<<field[i].x()<<' '<<field[i].y()<<' '<<field[i].z()<<tab;
+            }
+            else if constexpr( std::is_same_v<T,realx4>)
+            {
+                os << field[i].x() << ' ' << field[i].y() << ' ' << field[i].z() << ' ' << field[i].w() << tab;
+            }
+            else
+            {
+                os<<field[i]<<tab;
+            }
+            
+        }
+        else
+        {
+            if constexpr(std::is_same_v<T,realx3>)
+            {
+                os<<defValue.x()<<' '<<defValue.y()<<' '<<defValue.z()<<tab;
+            }
+            else if constexpr( std::is_same_v<T,realx4>)
+            {
+                os << defValue.x() << ' ' << defValue.y() << ' ' << defValue.z() << ' ' << defValue.w() << tab;
+            }
+            else
+            {
+                os<<defValue<<tab;
+            }
+        }
+
+    }
+    os<<endl;
+    return true;
+}
+
+}
+
+pFlow::postprocessOperation::postprocessOperation
+(
+    const dictionary &opDict, 
+    const regionPoints& regPoints,
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation
+    (
+        opDict,
+        opDict.getVal<word>("field"),
+        opDict.getValOrSet<word>("phi", "one"),
+        opDict.getValOrSet<word>("includeMask", "all"),
+        regPoints,
+        fieldsDB
+    )
+{}
+
+pFlow::postprocessOperation::postprocessOperation
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName,
+    const word& includeName,  
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    operationDict_(opDict),
+    threshold_
+    (
+        opDict.getValOrSet<int>("threshold", 1)
+    ),
+    divideByVolume_
+    (
+        opDict.getValOrSet<Logical>("divideByVolume", Logical(false))
+    ),
+    regionPoints_
+    (
+        regPoints
+    ),
+    database_
+    (
+        fieldsDB
+    ),
+    fieldName_
+    (
+        fieldName
+    ),
+    phiFieldName_
+    (
+        phiName
+    ),
+    includeMask_
+    (
+        includeMask::create(includeName, opDict, fieldsDB)
+    )
+{
+
+    if(!fieldsDB.getFieldType(fieldName_, fieldType_))
+    {
+        fatalErrorInFunction;
+        fatalExit;
+    }
+}
+const pFlow::Time& pFlow::postprocessOperation::time() const
+{
+    return database_.time();
+}
+
+bool pFlow::postprocessOperation::write(const fileSystem &parDir) const
+{   
+    auto ti = time().TimeInfo();
+
+    if(!osPtr_)
+    {
+        fileSystem path = parDir+(
+            processedFieldName() + ".Start_" + ti.timeName());
+        osPtr_ = makeUnique<oFstream>(path);
+        
+        regPoints().write(osPtr_());
+    }
+
+    const auto& field = processedField();
+
+    std::visit
+    (
+        [&](auto&& arg)->bool
+        {
+            return writeField(osPtr_(), ti.t(), arg, threshold_);
+        },
+        field
+    );
+
+    return true;
+}
+
+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/postprocessOperation.hpp

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

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -0,0 +1,162 @@
+#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()
+    ),
+    volumeFactor_
+    (
+        "volumeFactor",
+        regionPointsPtr_(),
+        1.0
+    ),
+    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( !forceUpdate && !timeControl().eventTime(ti))
+    {
+        executed_ = false;
+        return true;
+    }
+
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
+    // update processing methods
+    auto& regPoints = this->regPoints();
+    
+    if(!regPoints.update())
+    {
+        fatalErrorInFunction
+            << "regionPoints update failed  for "
+            << operationDicts_.globalName() 
+            << endl;
+        return false;
+    }
+
+    auto centers   = regPoints.centers();
+    const uint32 n = centers.size();
+    auto points    = this->database().updatePoints();
+    
+    for(uint32 i=0; i<n; i++)
+    {
+        auto indices = 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, volumeFactor_) )
+        {
+            fatalErrorInFunction
+                <<"error occured in executing operatoin defined in dict "
+                << op->operationDict()
+                <<endl;
+            return false;
+        }
+    }
+
+    executed_ = true;
+
+    return true;
+}
+
+template <typename RegionType, typename ProcessMethodType>
+inline 
+bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::write
+(
+    const fileSystem &parDir
+) const
+{
+    if(!executed_) return true;
+    
+    if(regionPointsPtr_().writeToSameTimeFile())
+    {
+        for(auto& operation:operatios_)
+        {
+            if(!operation->write(parDir))
+            {
+                fatalErrorInFunction
+                    <<"Error occurred in writing operation defined in dict "
+                    << operation->operationDict()
+                    <<endl;
+                return false;
+            }
+        }
+    }
+    else
+    {
+        notImplementedFunction;
+        return false;
+    }
+    
+
+
+    return true;
+}

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp

@@ -0,0 +1,131 @@
+/*------------------------------- 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_;
+
+    regionField<real>               volumeFactor_;
+
+    bool    executed_{false};
+
+    dictionaryList                  operationDicts_;
+
+protected:
+    
+    std::vector<ProcessMethodType>& regionProecessMethod()
+    {
+        return regionsProcessMethod_;
+    }
+
+    regionField<real>& volumeFactor()
+    {
+        return volumeFactor_;
+    }
+    
+    const regionField<real>& volumeFactor()const
+    {
+        return volumeFactor_;
+    }
+
+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_;
+    }
+
+    bool write(const fileSystem& parDir)const override;
+
+};
+
+
+
+}
+
+#include "PostprocessComponent.cpp"
+
+#endif

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp

@@ -0,0 +1,75 @@
+/*------------------------------- 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 __PostprocessComponentArithmetic_hpp__
+#define __PostprocessComponentArithmetic_hpp__
+
+#include "PostprocessComponent.hpp"
+#include "arithmetic.hpp"
+
+namespace pFlow
+{
+
+template<typename RegionType>
+class PostprocessComponentArithmetic
+:
+    public PostprocessComponent<RegionType, arithmetic>
+{
+public:
+
+    /// type info
+    TypeInfoTemplate12("PostprocessComponent", RegionType, arithmetic);
+
+    PostprocessComponentArithmetic  
+    (
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl
+    )
+    :
+        PostprocessComponent<RegionType, arithmetic>(dict, fieldsDB, defaultTimeControl)
+    {
+        /// initializes the arithmetic 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();
+        auto& volFactor = this->volumeFactor();
+        const uint32 n = d.size();
+        for(uint32 i=0; i<n; i++)
+        {
+            regs[i] = arithmetic(); // Changed from uniformDistribution() to arithmetic()
+            volFactor[i] = 1.0;
+        }
+    }
+
+    // add the virtual constructor
+    add_vCtor(
+        postprocessComponent,
+        PostprocessComponentArithmetic,
+        dictionary
+    );
+
+};
+
+}
+
+
+#endif //__PostprocessComponentArithmetic_hpp__

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp

@@ -0,0 +1,79 @@
+/*------------------------------- 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"
+#include "numericConstants.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();
+        auto& volFactor = this->volumeFactor();
+
+        const uint32 n = d.size();
+        for(uint32 i=0; i<n; i++)
+        {
+            auto r = d[i]/2;
+            regs[i] = GaussianDistribution(c[i], pow(r/3.0,2));
+            volFactor[i] = 0.677683 / (4.0/3.0*Pi*r);
+        }
+    }
+
+    // add the virtual constructor
+    add_vCtor
+    (
+        postprocessComponent,
+        PostprocessComponentGaussian,
+        dictionary
+    );
+
+};
+
+}
+
+
+#endif //__PostprocessComponentGaussian_hpp__

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp

@@ -0,0 +1,75 @@
+/*------------------------------- 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 __PostprocessComponentUniform_hpp__
+#define __PostprocessComponentUniform_hpp__
+
+#include "PostprocessComponent.hpp"
+#include "uniformDistribution.hpp"
+
+namespace pFlow
+{
+
+template<typename RegionType>
+class PostprocessComponentUniform
+:
+    public PostprocessComponent<RegionType, uniformDistribution>
+{
+public:
+
+    /// type info
+    TypeInfoTemplate12("PostprocessComponent", RegionType, uniformDistribution);
+
+    PostprocessComponentUniform
+    (
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl
+    )
+    :
+        PostprocessComponent<RegionType,uniformDistribution>(dict, fieldsDB, defaultTimeControl)
+    {
+        /// initializes the Uniform 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();
+        auto& volFactor = this->volumeFactor();
+        const uint32 n = d.size();
+        for(uint32 i=0; i<n; i++)
+        {
+            regs[i] = uniformDistribution();
+            volFactor[i] = 1.0;
+        }
+    }
+
+    // add the virtual constructor
+    add_vCtor(
+        postprocessComponent,
+        PostprocessComponentUniform,
+        dictionary
+    );
+
+};
+
+}
+
+
+#endif //__PostprocessComponentUniform_hpp__

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp

@@ -0,0 +1,42 @@
+/*------------------------------- 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"
+#include "PostprocessComponentUniform.hpp"
+#include "PostprocessComponentArithmetic.hpp"
+
+// region types 
+#include "sphereRegionPoints.hpp"
+#include "lineRegionPoints.hpp"
+#include "multipleSpheresRegionPoints.hpp"
+
+
+template class pFlow::PostprocessComponentGaussian<pFlow::sphereRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::sphereRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::sphereRegionPoints>;
+
+template class pFlow::PostprocessComponentGaussian<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::multipleSpheresRegionPoints>;
+
+template class pFlow::PostprocessComponentGaussian<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::lineRegionPoints>;
+

src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp

@@ -0,0 +1,152 @@
+#include "particleProbePostprocessComponent.hpp"
+#include "Time.hpp"
+
+namespace pFlow
+{
+    
+template<typename T>
+inline
+regionField<T> porbeExecute
+(
+    const word&          regFieldName,
+    const span<T>&       field,
+    const regionPoints&  regPoints
+)
+{
+    regionField<T> processedField(regFieldName, regPoints, T{});
+    auto partIndices = regPoints.indices(0);
+
+    uint n = 0;
+    for(auto index:partIndices)
+    {
+        if(index != -1)
+        {
+            processedField[n] =  field[index];
+        }
+        n++;
+    }
+    return processedField;
+}
+
+template<typename T>
+inline bool writeField
+(
+    iOstream&             os,
+    timeValue             t,
+    const regionField<T>& field,
+    const regionPoints&   regPoints,
+    const T&              invalidVal = T{}
+)
+{
+    auto indices = regPoints.indices(0);
+    const uint32 s= field.size();
+    os<< t <<tab;
+    for(uint32 i=0; i<s; i++)
+    {
+        if constexpr(std::is_same_v<T,realx3>)
+        {
+            os<<field[i].x()<<' '<<field[i].y()<<' '<<field[i].z()<<tab;
+        }
+        else if constexpr( std::is_same_v<T,realx4> )
+        {
+            os<<field[i].x()<<' '<<field[i].y()<<' '<<field[i].z()<<tab<<field[i].w()<<tab;
+        }
+        else
+        {
+            os<<field[i]<<tab;
+        }
+    }
+    os<<endl;
+    return true;
+}
+
+}
+
+pFlow::particleProbePostprocessComponent::particleProbePostprocessComponent
+(
+    const dictionary &dict, 
+    fieldsDataBase &fieldsDB, 
+    const baseTimeControl &defaultTimeControl
+)
+:
+    postprocessComponent(dict, fieldsDB, defaultTimeControl),
+    regionPointsPtr_
+    (
+        makeUnique<centerPointsRegionPoints>(dict, fieldsDB)
+    ),
+    fieldName_
+    (
+        dict.getVal<word>("field")
+    ),
+    name_(dict.name())
+{}
+
+bool pFlow::particleProbePostprocessComponent::execute
+(
+    const timeInfo &ti, 
+    bool forceExecute
+)
+{
+    if( !forceExecute && !timeControl().eventTime(ti))
+    {
+        executed_ = false;
+        return true;
+    }
+
+    REPORT(1)<<"Executing postprocess component ("
+        <<Blue_Text(ti.timeName())<<" s) : "
+        << name()
+        <<END_REPORT;
+
+    if(!regionPointsPtr_().update())
+    {
+        fatalErrorInFunction
+            << "regionPoints update for "<< name_ << " failed. \n";
+        return false;
+    }
+
+    auto field = database().updateFieldAll(fieldName_);
+    auto pFieldName = name_;
+    
+    processedField_ = makeUnique<processedRegFieldType>
+    (
+        std::visit(
+            [&](auto&& f) -> processedRegFieldType
+            {
+                return porbeExecute(
+                    pFieldName,
+                    f,
+                    regionPointsPtr_());
+            },
+            field)
+    );
+
+    executed_ = true;
+    return true;
+}
+
+
+bool pFlow::particleProbePostprocessComponent::write(const fileSystem& parDir)const
+{
+    if(! executed_ ) return true;
+    
+    const auto ti = database().time().TimeInfo();
+    
+    if( !osPtr_)
+    {
+        // file is not open yet
+        fileSystem path = parDir + (name_+".Start_"+ti.timeName());
+        osPtr_ = makeUnique<oFstream>(path);
+        regionPointsPtr_().write(osPtr_());
+    }
+
+    std::visit
+    ([&](auto&& arg)->bool
+        {
+            return writeField(osPtr_(), ti.t(), arg, regionPointsPtr_());
+        },
+        processedField_()
+    );
+
+    return true;
+}

src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.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 __particleProbePostprocessComponent_hpp__
+#define __particleProbePostprocessComponent_hpp__
+
+#include "postprocessComponent.hpp"
+#include "fieldsDataBase.hpp"
+#include "centerPointsRegionPoints.hpp"
+#include "regionField.hpp"
+#include "oFstream.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                                fieldName_; 
+
+    word                                name_;
+
+    mutable uniquePtr<oFstream>         osPtr_ = nullptr;
+
+public:
+
+    TypeInfo("PostprocessComponent<centerPoints,particleProbe>");
+
+    particleProbePostprocessComponent
+    (
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl
+    );
+    
+    ~particleProbePostprocessComponent()override = default;
+
+    add_vCtor
+    (
+        postprocessComponent,
+        particleProbePostprocessComponent,
+        dictionary
+    );
+     
+    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_;
+    }
+
+    bool write(const fileSystem& parDir)const override;
+
+};
+    
+}
+
+#endif //__particleProbePostprocessComponent_hpp__

src/PostprocessData/postprocessComponent/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/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 "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;
+
+    virtual 
+    bool write(const fileSystem& parDir)const = 0; 
+
+    static
+    uniquePtr<postprocessComponent> create(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDB,
+        const baseTimeControl& defaultTimeControl);
+
+};
+
+} // namespace pFlow
+
+#endif // __postprocessComponent_hpp__

src/PostprocessData/postprocessData/postprocessData.cpp

@@ -0,0 +1,137 @@
+
+/*------------------------------- 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 "postprocessGlobals.hpp"
+#include "postprocessComponent.hpp"
+
+pFlow::postprocessData::postprocessData(const systemControl &control)
+:
+    auxFunctions(control),
+    time_(control.time()),
+    fieldsDataBasePtr_
+    ( 
+        fieldsDataBase::create
+        (
+            const_cast<systemControl&>(control),
+            true
+        )
+    ),
+    dict_
+    (
+        objectFile
+        (
+            "postprocessDataDict",
+            control.settings().path(),
+            objectFile::READ_IF_PRESENT,
+            objectFile::WRITE_NEVER
+        )
+    )
+{
+    postProcessGlobals::defaultDir__ = CWD()/pFlow::postProcessGlobals::defaultRelDir__;
+    
+    // if dictionary is not provided, no extra action is required. 
+    if( !dict_.fileExist() || !dict_.headerOk() )
+    {
+        WARNING<<"You requested postprocessData function while,"
+            <<" the dictionary system/postprocessDataDict does not exist."
+            <<" This feature is disabled in the current run."<<END_WARNING;
+        return;
+    }
+
+    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");
+    }
+
+    shapeType_ = dict_.getValOrSet<word>
+    (
+        "shapeType", 
+        word("sphere")
+    );
+    
+    componentsDictsPtr_ = makeUnique<dictionaryList>(readDictList("components", dict_));
+
+    for(auto& compDict:*componentsDictsPtr_)
+    {
+        postprocesses_.push_back( postprocessComponent::create(
+            compDict, 
+            fieldsDataBasePtr_(), 
+            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;
+}
+
+bool pFlow::postprocessData::write() const
+{
+    for(auto& component:postprocesses_)
+    {   
+        if(!component->executed())
+        {
+            continue;
+        }
+        if(!component->write(postProcessGlobals::defaultDir__/component->name()))
+        {
+            fatalErrorInFunction
+                <<"Error occured in writing postprocess component: "
+                <<component->name()<<endl;
+            fatalExit;
+        }
+    }
+    return true;
+}

src/PostprocessData/postprocessData/postprocessData.hpp

@@ -0,0 +1,104 @@
+/*------------------------------- 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 
+    uniquePtr<fieldsDataBase>        fieldsDataBasePtr_;
+
+    /// file dictionary that is constructed from the file (postProcessDataDict)
+    fileDictionary        dict_;
+
+    /// name of the shape for use in the time of postprocess after simulation 
+    word                  shapeType_;
+
+    /// list of dictionaries for postprocess components 
+    uniquePtr<dictionaryList>        componentsDictsPtr_ = nullptr;
+
+    /// @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 write()const override;
+    
+
+};
+
+} // namespace pFlow
+
+#endif // __postprocessData_hpp__

src/PostprocessData/postprocessData/postprocessGlobals.hpp

@@ -0,0 +1,35 @@
+/*------------------------------- 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 __postProcessGlobals_hpp__
+#define __postProcessGlobals_hpp__
+
+#include "fileSystem.hpp"
+
+namespace pFlow::postProcessGlobals
+{
+
+static fileSystem defaultDir__;
+inline const word defaultRelDir__ = "postprocessData";
+
+}
+
+
+#endif // __postProcessGlobals_hpp__

src/PostprocessData/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/processMethod/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"
+#include "numericConstants.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_))/sqrt(2.0*Pi*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/processMethod/arithmetic.hpp

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

src/PostprocessData/processMethod/uniformDistribution.hpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __uniformDistribution_hpp__
+#define __uniformDistribution_hpp__
+
+#include <vector>
+
+#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<uint32>& indices,  
+        const span<realx3>& points
+    )
+    {
+        return updateWeights(indices.size());
+    }
+
+    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/region/regionFields/regionField.hpp

@@ -0,0 +1,126 @@
+/*------------------------------- 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
+{
+public:
+    using FieldType         = Field<T, HostSpace>;
+
+    using iterator        	= typename FieldType::iterator;
+
+    using const_iterator   	= typename FieldType::const_iterator;
+
+    using reference       	= typename FieldType::reference;
+    
+    using const_reference  	= typename FieldType::const_reference;
+
+    using value_type       	= typename FieldType::value_type;
+    
+    using pointer         	= typename FieldType::pointer;
+    
+    using const_pointer    	= typename FieldType::const_pointer;
+private:
+           
+    /// the field value
+    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;
+
+    const regionPoints& regPoints() const
+    {
+        return regionPoints_;
+    }
+
+    /// 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();
+    }
+
+    auto size()const
+    {
+        return field_.size();
+    }
+
+    bool empty()const
+    {
+        return field_.empty();
+    }
+
+    /// @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/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/centerPointsRegionPoints.cpp

@@ -0,0 +1,97 @@
+#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
+    // TODO: this should be corrected to select ids of particles
+    // that are selected based on the selector (this is visa versa)
+    {
+        auto selectorPtr = pStructSelector::create(
+            selector, 
+            database().pStruct(), 
+            probDict_.subDict(selector+"Info"));
+        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 true;
+}
+
+pFlow::centerPointsRegionPoints::centerPointsRegionPoints(
+    const dictionary &dict,
+    fieldsDataBase &fieldsDataBase)
+    : regionPoints(dict, fieldsDataBase),
+      idName_(dict.getValOrSet<word>("idName", "id")),
+      probDict_(dict)
+{}
+
+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 true;
+}
+
+bool pFlow::centerPointsRegionPoints::write(iOstream &os) const
+{
+    if(firstTimeUpdate_)
+    {
+        fatalErrorInFunction
+            <<"Write operation should be done before executing the update. \n";
+        return false;   
+    }
+    os <<"# Probing particles for selected ids\n"; 
+    
+    os << "# SelectedId: ";
+
+    for(auto& id: ids_)
+    {
+        os << id << tab;
+    }
+    os << endl;
+    os << "time"<<endl;
+
+    return true;
+}

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

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

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

@@ -0,0 +1,120 @@
+#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);
+    diameters_.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());
+}
+
+pFlow::span<pFlow::uint32> pFlow::lineRegionPoints::indices(uint32 elem)
+{
+    if(elem >= size())
+    {
+        fatalErrorInFunction
+            << "The element index is out of range. elem: " << elem
+            << " size: " << size() << endl;
+        fatalExit;
+    }
+    
+    return span<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;
+}
+
+bool pFlow::lineRegionPoints::write(iOstream &os) const
+{
+    os << "# Spheres along a straight line \n";
+    os << "#  No." << tab << "centerPoint" << tab << "diameter" << endl;
+    for(uint32 i=0; i < sphereRegions_.size(); ++i)
+    {
+        os << "# " << i << tab << sphereRegions_[i].center() 
+           << tab << diameters_[i] << '\n';
+    }
+    os << "time/No. ";
+    for(uint32 i=0; i < sphereRegions_.size(); ++i)
+    {
+        os << i << tab;
+    }
+    os << endl;
+    return true;
+}

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

@@ -0,0 +1,148 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+/**
+ * @class lineRegionPoints 
+ * @brief Spherical regions along a line for selecting points/particles
+ * 
+ * The lineRegionPoints class is responsible for selecting points/particles along a 
+ * specified line and creating sphere regions around those points. It partitions
+ * the line into multiple sphere regions (equally spaced) and keeps track of 
+ * which points/particles in the simulation fall into each region.
+ * 
+ * This class is used for post-processing data by analyzing distributions of
+ * particles along a linear path through the simulation domain. It maintains:
+ *   - A line defining the sampling path
+ *   - Spherical regions along this line that include particles
+ *   - Center points for each region
+ *   - Volumes and diameters of regions
+ *   - Indices of points/particles contained in each region
+ * 
+ * The regions can be updated as the simulation progresses, and the data
+ * can be written to output for analysis.
+ *
+ * @see regionPoints
+ * @see line
+ * @see sphere
+ * @see fieldsDataBase
+ */
+
+#ifndef __lineRegionPoints_hpp__
+#define __lineRegionPoints_hpp__
+
+
+#include "regionPoints.hpp"
+#include "line.hpp"
+#include "sphere.hpp"
+#include "Vectors.hpp"
+
+
+namespace pFlow
+{
+
+class lineRegionPoints
+:
+    public regionPoints
+{
+private:
+
+    /// Line path defining the axis of the spherical regions
+    line                line_;
+
+    /// Collection of sphere regions along the line
+    Vector<sphere>        sphereRegions_;
+
+    /// Center points of all spherical regions
+    realx3Vector          centerPoints_;
+
+    /// Volumes of all spherical regions
+    realVector            volumes_;
+
+    /// Diameters of all spherical regions
+    realVector            diameters_;
+
+    /// Point/particles indices selected by each region
+    Vector<uint32Vector>  selectedPoints_;
+
+public:
+
+    /// Type information for runtime type identification
+    TypeInfo(line::TYPENAME());
+
+    /// Construct from dictionary that contains lineInfo and fields database
+    lineRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    /// Default destructor
+    ~lineRegionPoints() override = default;
+
+    /// Return number of regions
+    uint32 size()const override
+    {
+        return sphereRegions_.size();
+    }
+
+    /// Check if regions list is empty
+    bool empty()const override
+    {
+        return sphereRegions_.empty();
+    }
+
+    /// Return volumes of all regions
+    span<const real> volumes()const override
+    {
+        return span<const real>(volumes_.data(), volumes_.size());
+    }
+
+    /// Return equivalent diameters of all regions
+    span<const real> eqDiameters()const override
+    {
+        return span<const real>(diameters_.data(), diameters_.size());
+    }
+
+    /// Return center points of all regions
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(centerPoints_.data(), centerPoints_.size());
+    }
+
+    /// Return indices of points in the specified element/region
+    span<const uint32> indices(uint32 elem)const override;
+
+    span<uint32> indices(uint32 elem) override;
+    
+    /// Update regions based on current particle positions
+    bool update() override;
+
+    /// Whether to write all data to the same time file
+    bool writeToSameTimeFile()const override
+    {
+        return true;
+    }
+
+    /// Write data to output stream
+    bool write(iOstream& os) const override;
+};
+    
+}
+
+
+#endif // __lineRegionPoints_hpp__

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

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

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

@@ -0,0 +1,166 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+/**
+ * @class       multipleSpheresRegionPoints
+ * @brief       A class to select and track particles contained within multiple 
+ *              spherical regions
+ * @details     This class defines multiple spherical regions in the simulation 
+ *             domain and identifies which particles are contained within each 
+ *             sphere at each time step. It inherits from the regionPoints base 
+ *             class and specializes it for handling multiple spherical regions 
+ *             simultaneously.
+ *             
+ *             The class reads a list of center points and radii from a dictionary,
+ *             creates sphere objects for each, and provides methods to:
+ *             - Track which particles are inside each spherical region
+ *             - Return volumetric information about the regions
+ *             - Access center points and size information for each sphere
+ *
+ *             
+ * @note        Used in post-processing workflows to analyze particle behavior
+ *             in specific regions of interest within the simulation domain.
+ * 
+ * @see        regionPoints Base class for all region-based point selection
+ * @see        sphere Geometric primitive used to define spherical regions
+ * @see        postprocessPhasicFlow Utility for post-processing simulation data
+ * @see        fieldsDataBase Database containing simulation field data
+ */
+
+#ifndef __multipleSpheresRegionPoints_hpp__
+#define __multipleSpheresRegionPoints_hpp__
+
+#include "regionPoints.hpp"
+#include "sphere.hpp"
+#include "Vectors.hpp"
+
+namespace pFlow
+{
+
+class multipleSpheresRegionPoints
+:
+    public regionPoints
+{
+private:
+
+    /// Vector containing all spherical regions used for particle selection
+    Vector<sphere>        sphereRegions_;
+
+    /// Center coordinates of all spherical regions
+    realx3Vector          centerPoints_;
+
+    /// Diameters of all spherical regions
+    realVector            diameters_;
+
+    /// Volumes of all spherical regions
+    realVector            volumes_;
+
+    /// Vectors of point indices for particles contained in each spherical region
+    /// Each element corresponds to a particular sphere region
+    Vector<uint32Vector>  selectedPoints_;
+
+public:
+
+    /// Type identification for run-time type information
+    TypeInfo("multipleSpheres");
+
+    /// Constructor
+    /// @param dict Dictionary containing multipleSpheresInfo for the regions
+    /// @param fieldsDataBase Reference to the database containing field data
+    multipleSpheresRegionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+
+    /// Virtual destructor for proper inheritance cleanup
+    ~multipleSpheresRegionPoints() override = default;
+
+    /// Returns the number of spherical regions
+    /// @return Number of spherical regions
+    uint32 size()const override
+    {
+        return sphereRegions_.size();
+    }
+
+    /// Checks if there are any spherical regions defined
+    /// @return True if no regions exist, false otherwise
+    bool empty()const override
+    {
+        return sphereRegions_.empty();
+    }
+
+    /// Returns the volumes of all spherical regions
+    /// @return Span containing the volumes of all regions
+    span<const real> volumes()const override
+    {
+        return span<const real>(volumes_.data(), volumes_.size());
+    }
+
+    /// Returns the diameters of all spherical regions
+    /// @return Span containing the diameters of all regions
+    span<const real> diameters()const
+    {
+        return span<const real>(diameters_.data(), diameters_.size());
+    }
+
+    /// Returns the equivalent diameters of all spherical regions
+    /// @return Span containing the equivalent diameters (same as diameters)
+    span<const real> eqDiameters()const
+    {
+        return diameters();
+    }
+
+    /// Returns the center coordinates of all spherical regions
+    /// @return Span containing the center points of all regions
+    span<const realx3> centers()const override
+    {
+        return span<const realx3>(centerPoints_.data(), centerPoints_.size());
+    }
+
+    /// Returns the indices of particles contained in a specific spherical region
+    /// @param elem Index of the spherical region to query
+    /// @return Span containing indices of particles within the specified region
+    span<const uint32> indices(uint32 elem)const override;
+
+    /// Returns the indices of particles contained in a specific spherical region
+    /// @param elem Index of the spherical region to query
+    /// @return Span containing indices of particles within the specified region
+    span<uint32> indices(uint32 elem) override;
+
+    /// Updates the selection of particles within each spherical region
+    /// @return True if update was successful, false otherwise
+    bool update() override;
+
+    /// Determines if data should be written to the same time file
+    /// @return True to indicate regions should be written to the same time file
+    bool writeToSameTimeFile()const override
+    {
+        return true;
+    }
+
+    /// Writes region data to the output stream
+    /// @param os Output stream to write data to
+    /// @return True if write operation was successful, false otherwise
+    bool write(iOstream& os) const override;
+};
+    
+}
+
+
+#endif // __multipleSpheresRegionPoints_hpp__

src/PostprocessData/region/regionPoints/regionPoints/regionPoints.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/regionPoints.hpp

@@ -0,0 +1,118 @@
+/*------------------------------- 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 "iOstream.hpp"
+#include "dictionary.hpp"
+#include "pointStructure.hpp"
+
+
+namespace pFlow
+{
+
+class fieldsDataBase;
+class Time;
+
+/**
+ * @class regionPoints
+ * @brief Abstract base class for managing and processing volumetric regions 
+ * in the simulation. Particles are selected based on their positions within
+ * these defined regions.
+ * 
+ * This class provides an interface for accessing and manipulating data 
+ * related to regions of points (particles), including their volumes, equivalent 
+ * diameters, center points, and particle indices that they contain. It interacts with the 
+ * fieldsDataBase and Time classes to retrieve simulation data. 
+ */
+class regionPoints
+{
+        using PointsTypeHost = typename pointStructure::PointsTypeHost;
+
+    /// Reference to the fields database containing simulation data
+    fieldsDataBase& fieldsDataBase_;
+
+public:
+
+    TypeInfo("regionPoints");
+
+    /// Constructor with dictionary and fields database reference
+    regionPoints(
+        const dictionary& dict,
+        fieldsDataBase& fieldsDataBase);
+    
+    /// Default virtual destructor
+    virtual ~regionPoints() = default;
+
+    /// Returns reference to the time object from the database
+    const Time& time()const;
+
+    /// Returns const reference to the fields database
+    const fieldsDataBase& database()const;
+
+    /// Returns non-const reference to the fields database
+    fieldsDataBase& database();
+    
+    /// @brief size of elements 
+    virtual 
+    uint32 size()const = 0;
+    
+    /// @brief  check if the region is empty
+    virtual 
+    bool empty()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;
+
+    /// Returns const span of particle indices inside the specified element region
+    virtual 
+    span<const uint32> indices(uint32 elem)const = 0;
+
+    /// Returns non-const span of particle indices inside the specified element region
+    virtual
+    span<uint32> indices(uint32 elem) = 0;
+
+    /// Updates the points (particles) inside regions based on current particle positions
+    virtual 
+    bool update() = 0;
+    
+    /// Returns true if the region should be written to the same time file
+    virtual
+    bool writeToSameTimeFile()const = 0;
+
+    /// Writes region data to the output stream
+    virtual 
+    bool write(iOstream& os)const=0;
+    
+};
+
+}
+
+#endif // __regionPoints_hpp__

src/PostprocessData/region/regionPoints/sphereRegionPoints/sphereRegionPoints.cpp

@@ -0,0 +1,41 @@
+#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 true;
+}
+
+bool pFlow::sphereRegionPoints::write(iOstream &os) const
+{
+    os <<"# Single sphere\n";
+    os <<"# center point: "<<sphereRegion_.center()<<endl;
+    os <<"# diameter: "<< diameter_ << endl;
+    os <<"time"<< tab <<"value"<<endl;
+    
+    return true;
+}

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

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

src/PostprocessData/sampleDictionary/postprocessDataDict

@@ -0,0 +1,129 @@
+/* -------------------------------*- C++ -*--------------------------------- *\ 
+|  phasicFlow File                                                            | 
+|  copyright: www.cemf.ir                                                     | 
+\* ------------------------------------------------------------------------- */  
+objectName   postprocessDataDict;
+objectType   dictionary;;
+fileFormat   ASCII;
+/*---------------------------------------------------------------------------*/
+
+runTimeActive yes;
+
+defaultTimeControl
+{
+	timeControl         timeStep;
+	startTime           0;
+	endTime             1000;
+	executionInterval   150; 
+}
+
+components
+(
+      velocityProb
+      {
+          processMethod    particleProbe;
+          processRegion    centerPoints;
+          selector        id;
+          field   component(position,y);
+          ids (0 10 100);
+      }
+    
+      onSingleSphere
+      {
+          // method of performing the sum (arithmetic, uniformDistribution, GaussianDistribution)
+          processMethod    arithmetic; 
+          processRegion    sphere; // type of region on which processing is performed 
+          
+          sphereInfo
+          {
+              radius 0.01;
+              center (-0.08 -0.08 0.015);
+          }
+          
+          timeControl default; // settings, timeStep, simulationTime 
+          
+          /// all the post process operations to be done 
+          operations
+          (
+              // computes the arithmetic mean of particle velocity 
+              averageVel
+              {
+                  function	  average;
+                  field           velocity;
+                  divideByVolume 	  no;  //default 
+                  threshold 	  3;   //default is 1;
+                  includeMask     all;
+              }
+              
+              // computes the fraction of par1 in the region 
+              par1Fraction
+              {
+                  function        average;
+                  field           one;
+                  phi             one; // default
+                  divideByVolume  no;
+                  includeMask     lessThan;
+                  
+                  // diameter of par1 is 0.003, so these settings 
+                  // will select only particles of type par1
+                  lessThanInfo
+                  {
+                      field      diameter;
+                      
+                      value      0.0031;  
+                  }
+              }
+
+              numberDensity
+              {
+                  function        sum;
+                  field           one;
+                  phi             one; // default
+                  divideByVolume  yes;
+                  
+              }
+          );
+    }
+    
+    alongALine
+    {
+          processMethod    arithmetic; 
+          processRegion    line;
+          
+          // the time interval for executing the post-processing 
+          // other options: timeStep, default, and settings
+          timeControl simulationTime; 
+          startTime   1.0;
+          endTime     3.0;
+          executionInterval 0.1;
+           
+          
+          // 10 spheres with radius 0.01 along the straight line defined by p1 and  p2 
+          lineInfo
+          {
+            p1 (0 0 0);
+            p2 (0 0.15 0.15);
+            numPoints 10;
+            radius    0.01;
+          }
+          
+          operations
+          (
+              // computes the arithmetic mean of particle velocity 
+              numberDensity
+              {
+                  function	   sum;
+                  field            one;
+                  divideByVolume  yes;  //default is no 
+              }
+              
+              volumeDensity
+              {
+                function         sum;
+                field            cube(diameter); // d^3, although it differs by pi/6
+                divideByVolume   yes;  //default is no 
+              }
+         );
+    }
+    
+);

src/phasicFlow/CMakeLists.txt

@@ -47,6 +47,7 @@ repository/Time/Time.cpp
 repository/Time/timeControl.cpp
 repository/Time/baseTimeControl.cpp
 repository/systemControl/systemControl.cpp
+repository/systemControl/auxFunctions/auxFunctions.cpp
 repository/systemControl/dynamicLinkLibs.cpp
 
 commandLine/commandLine.cpp

src/phasicFlow/containers/List/anyList/anyList.hpp

@@ -79,6 +79,27 @@ public:
 		/// List of varibales names
 		const wordList& names()const;
 
+		template<typename T, typename... Args>
+		void emplaceBackOrReplace(const word& name, Args&&... args)
+		{
+			if( contains(name))
+			{
+				int32 i = names_.findi(name);
+				types_[i] = getTypeName<T>();
+				anyList_.pos(i)->reset();
+				anyList_.pos(i)-> emplace<T>(std::forward<Args>(args)...);
+			}
+			else
+			{
+				names_.push_back(name);
+				types_.push_back(getTypeName<T>());
+				anyList_.emplace_back( 
+					std::in_place_type<T>, 
+					std::forward<Args>(args)...);
+			}
+			
+		}
+
 		/// Create variable using constructor in-place 
 		template<typename T, typename... Args>
 		reference emplaceBack(const word& name, Args&&... args)

src/phasicFlow/dictionary/dictionaryList.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 __dictionaryList_hpp__
+#define __dictionaryList_hpp__
+
+#include "List.hpp"
+#include "dictionary.hpp"
+
+
+namespace pFlow
+{
+
+class dictionaryList   
+:
+    public List<dictionary>
+{
+private:
+
+    const dictionary& parrent_;
+
+    word              name_;
+
+
+public:
+    
+    TypeInfoNV("dictionaryList");
+
+    dictionaryList(const word& name,  const dictionary& parrent)
+    :
+        List<dictionary>(),
+        parrent_(parrent),
+        name_(name)
+    {}
+    
+    dictionaryList(const dictionaryList&) = default;
+
+    dictionaryList(dictionaryList&&) = default;
+
+    dictionaryList& operator = (const dictionaryList&) = default;
+
+    dictionaryList& operator = (dictionaryList&&) = default;
+
+    ~dictionaryList() = default;
+
+    const word& name()const
+    {
+        return name_;
+    }
+
+    const word globalName()const
+    {
+        return parrent_.globalName();
+    }
+
+    List<dictionary>& dicts()
+    {
+        return *this;
+    }
+
+    const List<dictionary>& dicts()const
+    {
+        return *this;
+    }
+
+    const dictionary& parrent()const
+    {
+        return parrent_;
+    }
+    
+};
+
+inline
+dictionaryList readDictList(const word& name, const dictionary& dict)
+{
+    if(!dict.containsDataEntry(name))
+    {
+        fatalErrorInFunction
+            <<"data entry: "<< name
+            <<" does not exist in dictionary "<< dict.globalName()<<endl;
+    }
+
+    dictionaryList allDicts(name, dict);
+    auto dicts = dict.getVal<List<dictionary>>(name);
+    for(auto& d:dicts)
+    {
+        allDicts.emplace_back(d.name(), dict, d);
+    }
+
+    return allDicts;
+}
+
+}
+
+
+
+#endif

src/phasicFlow/globals/vocabs.hpp

@@ -47,6 +47,7 @@ const inline char* const contactSearchFile__  = "contactSearch";
 const inline char* const propertyFile__       = "interaction";
 const inline char* const interactionFile__    = "interaction";
 const inline char* const postprocessFile__    = "postprocessDict";
+const inline char* const postprocessDataFiel__ = "postprocessDataDict"; 
 
 const inline char* const uniform__    = "uniform";
 const inline char* const nonUniform__ = "nonUniform";

src/phasicFlow/repository/Time/baseTimeControl.cpp

@@ -21,14 +21,24 @@ Licence:
 #include "baseTimeControl.hpp"
 #include "timeInfo.hpp"
 
-pFlow::baseTimeControl::baseTimeControl
+void pFlow::baseTimeControl::setTimeControl
 (
-    const dictionary &dict,
-    const word& intervalPrefix,
-	timeValue defStartTime
+	timeValue startTime, 
+	timeValue endTime, 
+	timeValue interval, 
+	const word &intervalPrefix
 )
-:
-	intervalPrefix_(intervalPrefix)
+{
+	isTimeStep_ = false;
+	intervalPrefix_ = intervalPrefix;
+	rRange_ = stridedRange<timeValue>(startTime, endTime, interval);
+}
+
+pFlow::baseTimeControl::baseTimeControl(
+    const dictionary &dict,
+    const word &intervalPrefix,
+    timeValue defStartTime)
+    : intervalPrefix_(intervalPrefix)
 {
 	auto tControl = dict.getVal<word>("timeControl");
 	if(tControl == "timeStep")
@@ -119,6 +129,24 @@ pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, co
 {
 }
 
+pFlow::baseTimeControl::baseTimeControl
+(
+	timeValue start, 
+	timeValue end, 
+	timeValue stride, 
+	const word &intervalPrefix
+)
+:
+	isTimeStep_(false),
+	rRange_(
+		start, end, stride
+	),
+	intervalPrefix_(
+		intervalPrefix.size()==0uL? word("interval"): intervalPrefix+"Interval"
+	)
+{ 
+}
+
 bool pFlow::baseTimeControl::eventTime(uint32 iter, timeValue t, timeValue dt) const
 {
 	if(isTimeStep_)

src/phasicFlow/repository/Time/baseTimeControl.hpp

@@ -32,13 +32,21 @@ class baseTimeControl
 {
 private:
 
-	bool              isTimeStep_;
+    bool                     isTimeStep_;
 
-	int32StridedRagne iRange_;
+    int32StridedRagne        iRange_;
 
-	stridedRange<timeValue>  rRange_;
+    stridedRange<timeValue>  rRange_;
 
-	const word        intervalPrefix_;
+    word                      intervalPrefix_;
+
+protected:
+
+	void setTimeControl(
+		timeValue startTime, 
+		timeValue endTime, 
+		timeValue interval, 
+		const word& intervalPrefix);
 
 public:
 
@@ -59,8 +67,33 @@ public:
 		int32 end,
 		int32 stride,
 		const word& intervalPrefix = ""
+	); 
+
+	baseTimeControl(
+		timeValue start,
+		timeValue end,
+		timeValue stride,
+		const word& intervalPrefix = ""
 	);
 
+    baseTimeControl(
+        const baseTimeControl& other
+    ) = default;
+
+    baseTimeControl(
+        baseTimeControl&& other
+    ) = default;
+
+    baseTimeControl& operator=(
+        const baseTimeControl& other
+    ) = default;
+    
+    baseTimeControl& operator=(
+        baseTimeControl&& other
+    ) = default;
+    
+    ~baseTimeControl() = default;
+
 	inline bool isTimeStep() const
 	{
 		return isTimeStep_;

src/phasicFlow/repository/Time/timeControl.hpp

@@ -120,6 +120,16 @@ public:
 		return startTime_;
 	}
 
+	timeValue endTime()const
+	{
+		return endTime_;
+	}
+
+	timeValue saveInterval()const
+	{
+		return saveInterval_;
+	}
+
 	word timeName()const;	
 
 	timeValue currentTime() const 

src/phasicFlow/repository/Time/timeInfo.hpp

@@ -106,6 +106,12 @@ public:
 		return real2FixedStripZeros(currentTime_, presicion_);
 	}
 
+	inline 
+	word prevTimeName()const
+	{
+		return real2FixedStripZeros( max(currentTime_-dt_, timeValue(0)), presicion_);
+	}
+
 	static
 	uint32 precision()
 	{

src/phasicFlow/repository/repository/repository.cpp

@@ -236,7 +236,42 @@ size_t pFlow::repository::numRepositories()const
 	return repositories_.size();
 }
 
-pFlow::repository& pFlow::repository::lookupRepository(const word& name)
+const pFlow::IOobject *pFlow::repository::lookupObjectPtr
+(
+	const word &name
+) const
+{
+	if( auto [iter, success] = objects_.findIf(name); success )
+	{
+		return iter->second;
+	}
+	else
+	{
+		fatalErrorInFunction << 
+		"Object with name " << name << " is not found in repository " << this->name()<<endl <<
+		"list of avaiable objest is \n" << objectNames();
+		fatalExit;
+	}
+    return nullptr;
+}
+
+pFlow::IOobject *pFlow::repository::lookupObjectPtr(const word &name)
+{
+	if( auto [iter, success] = objects_.findIf(name); success )
+	{
+		return iter->second;
+	}
+	else
+	{
+		fatalErrorInFunction << 
+		"Object with name " << name << " is not found in repository " << this->name()<<endl <<
+		"list of avaiable objest is \n" << objectNames();
+		fatalExit;
+	}
+    return nullptr;
+}
+
+pFlow::repository &pFlow::repository::lookupRepository(const word &name)
 {
 
 	if( auto [iter, success] = repositories_.findIf(name); success )

src/phasicFlow/repository/repository/repository.hpp

@@ -168,6 +168,10 @@ public:
 			}
 		}
 		
+		const IOobject* lookupObjectPtr(const word& name)const;
+
+		IOobject* lookupObjectPtr(const word& name);
+
 		/// return a ref to the underlaying data in the object 
 		template<typename T>
 		T& lookupObject(const word& name);

src/phasicFlow/repository/systemControl/auxFunctions/auxFunctions.cpp

@@ -0,0 +1,56 @@
+/*------------------------------- 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 "auxFunctions.hpp"
+#include "systemControl.hpp"
+
+
+pFlow::uniquePtr<pFlow::auxFunctions> 
+    pFlow::auxFunctions::create(const systemControl& control)
+{
+    const auto& setDict = control.settingsDict();
+
+    if( !setDict.containsDataEntry("auxFunctions"))
+    {
+        return nullptr;
+    }
+
+    word func = setDict.getVal<word>("auxFunctions");
+
+    if( systemControlvCtorSelector_.search(func) )
+    {
+        REPORT(1)<<"Creating auxiliary function "<< Green_Text(func)<< " ..."<<END_REPORT;
+        return systemControlvCtorSelector_[func](control);
+    }
+    else
+    {
+        printKeys
+		( 
+			fatalError 
+                << "Ctor Selector "<< Yellow_Text(func) << " dose not exist "
+                << "for axuFunctions.\n"
+			    <<"Avaiable ones are: \n"
+			,
+			systemControlvCtorSelector_
+		);
+		fatalExit;
+        return nullptr;
+    }
+}

src/phasicFlow/repository/systemControl/auxFunctions/auxFunctions.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 __auxFunctions_hpp__
+#define __auxFunctions_hpp__
+
+#include "timeInfo.hpp"
+#include "uniquePtr.hpp"
+#include "virtualConstructor.hpp"
+
+namespace pFlow
+{
+
+class systemControl;
+
+class auxFunctions
+{
+
+public:
+
+    TypeInfo("auxFunctions");
+
+    auxFunctions(const systemControl& control){};
+
+    virtual ~auxFunctions()=default;
+
+    create_vCtor
+    (
+        auxFunctions,
+        systemControl,
+        (const systemControl& control),
+        (control)
+    );
+
+    virtual 
+    bool execute() = 0;
+    
+    virtual 
+    bool write()const = 0 ; 
+
+    static
+    uniquePtr<auxFunctions> create(const systemControl& control);
+
+};
+
+} // namespace pFlow
+
+#endif 

src/phasicFlow/repository/systemControl/systemControl.cpp

@@ -24,6 +24,10 @@ Licence:
 #include "iOstream.hpp"
 #include "types.hpp"
 #include "vocabs.hpp"
+#include "Logical.hpp"
+
+
+inline static bool axuFunctionsInitialized__ = false;
 
 bool pFlow::systemControl::readIncludeExclue(const dictionary& dict)
 {
@@ -185,7 +189,20 @@ pFlow::systemControl::systemControl(
 
 bool pFlow::systemControl::operator++(int)
 {
-	auto toContinue = time()++;
+	
+    if(!axuFunctionsInitialized__)
+    {
+        auxFunctions_ = auxFunctions::create(*this);
+        axuFunctionsInitialized__ = true;
+    }
+    
+    if(auxFunctions_)
+    {
+        auxFunctions_().execute();
+        auxFunctions_().write();
+    }
+
+    auto toContinue = time()++;
 
 	if (toContinue)
 	{
@@ -221,3 +238,11 @@ bool pFlow::systemControl::operator++(int)
 
 	return toContinue;
 }
+
+bool pFlow::systemControl::keepIntegrationHistory()const
+{
+    auto keepHistory = settingsDict_().getValOrSet(
+        "integrationHistory",
+        Logical{false});
+    return keepHistory();
+}

src/phasicFlow/repository/systemControl/systemControl.hpp

@@ -34,10 +34,13 @@ Licence:
 #include "Timers.hpp"
 #include "dynamicLinkLibs.hpp"
 #include "Set.hpp"
+#include "auxFunctions.hpp"
+
 
 namespace pFlow
 {
 
+
 class systemControl
 :
 	public repository
@@ -83,6 +86,7 @@ protected:
 
 	wordSet 		excludeList_;
     
+	uniquePtr<auxFunctions> 		auxFunctions_ = nullptr;
 	
 	bool readIncludeExclue(const dictionary& dict);
 
@@ -189,6 +193,7 @@ public:
 		return outFilePrecision_;
 	}
 
+	bool keepIntegrationHistory()const;
 	
 	bool isIncluded(const word& objName)const final
 	{

src/phasicFlow/smartPointers/uniquePtr.hpp

@@ -38,15 +38,16 @@ namespace pFlow
 
 
 template<
-	typename T
+	typename T,
+	typename Deleter = std::default_delete<T>
 >
 class uniquePtr
 :
-	public std::unique_ptr<T>
+	public std::unique_ptr<T, Deleter>
 {
 public:
 
-	using uniquePtrType = std::unique_ptr<T>;
+	using uniquePtrType = std::unique_ptr<T, Deleter>;
 
 	// using base constructors 
 	using uniquePtrType::unique_ptr;

src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp

@@ -93,9 +93,18 @@ pFlow::selectorRandomPoints::selectorRandomPoints(
   const dictionary&     dict
 )
   : pStructSelector(type, pStruct, dict),
-    begin_(dict.getVal<uint32>("begin")),
-    end_(dict.getValOrSet("end", pStruct.size())),
-    number_(dict.getValOrSet("number", 1))
+    begin_
+	(
+		dict.getValOrSet("begin",0u)
+	),
+    end_
+	(
+		dict.getValOrSetMax<uint32>("end", pStruct.size())
+	),
+    number_
+	(
+		dict.getValOrSetMax<uint32>("number", 0u)
+	)
 {
 	begin_  = max(begin_, 1u);
 	end_    = min(end_, static_cast<uint32>(pStruct.size()));

src/phasicFlow/structuredData/sphere/sphere.hpp

@@ -105,6 +105,12 @@ public:
 			return sqrt(radius2_);
 		}
 
+		INLINE_FUNCTION_HD
+		real volume()const
+		{
+			return 4.0/3.0* Pi * pow(radius(),3.0); 
+		}
+
 	//// - IO operation 
 		FUNCTION_H
 		bool read(iIstream & is);

src/phasicFlow/types/quadruple/quadruple.hpp

@@ -278,6 +278,6 @@ struct quadruple
 } // pFlow
 
 #include "quadrupleI.hpp"
-// #include "quadrupleMath.hpp"
+#include "quadrupleMath.hpp"
 
 #endif