diff --git a/src/Particles/GrainParticles/grainShape/grainShape.cpp b/src/Particles/GrainParticles/grainShape/grainShape.cpp
index 4558940b..246b5be8 100644
--- a/src/Particles/GrainParticles/grainShape/grainShape.cpp
+++ b/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_;
 }
diff --git a/src/Particles/GrainParticles/grainShape/grainShape.hpp b/src/Particles/GrainParticles/grainShape/grainShape.hpp
index e22a729d..6f0ad1f9 100644
--- a/src/Particles/GrainParticles/grainShape/grainShape.hpp
+++ b/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;
 
diff --git a/src/Particles/SphereParticles/sphereShape/sphereShape.cpp b/src/Particles/SphereParticles/sphereShape/sphereShape.cpp
index fac52d3c..8a5b478a 100644
--- a/src/Particles/SphereParticles/sphereShape/sphereShape.cpp
+++ b/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) )
diff --git a/src/Particles/SphereParticles/sphereShape/sphereShape.hpp b/src/Particles/SphereParticles/sphereShape/sphereShape.hpp
index b59625ce..10f308a1 100644
--- a/src/Particles/SphereParticles/sphereShape/sphereShape.hpp
+++ b/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;
diff --git a/src/Particles/particles/particles.hpp b/src/Particles/particles/particles.hpp
index 8c5f7741..b098f167 100644
--- a/src/Particles/particles/particles.hpp
+++ b/src/Particles/particles/particles.hpp
@@ -186,7 +186,7 @@ public:
 	}
 
 	inline 
-	uint maxId()const
+	uint32 maxId()const
 	{
 		return idHandler_().maxId();
 	}
diff --git a/src/Particles/particles/shape/shape.cpp b/src/Particles/particles/shape/shape.cpp
index bcf66e74..df84ff90 100644
--- a/src/Particles/particles/shape/shape.cpp
+++ b/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;
-}
\ No newline at end of file
+}
+
+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;
+	}
+
+}
diff --git a/src/Particles/particles/shape/shape.hpp b/src/Particles/particles/shape/shape.hpp
index 5a4007aa..4c7d74e9 100644
--- a/src/Particles/particles/shape/shape.hpp
+++ b/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);
+
 };
 
 }
diff --git a/src/PostprocessData/CMakeLists.txt b/src/PostprocessData/CMakeLists.txt
index bb3999d4..a68a6d75 100644
--- a/src/PostprocessData/CMakeLists.txt
+++ b/src/PostprocessData/CMakeLists.txt
@@ -4,6 +4,7 @@ set(SourceFiles
 
     # Fields database
     fieldsDataBase/fieldsDataBase.cpp
+    fieldsDataBase/simulationFieldsDataBase.cpp
   
     # Regions
     region/regionPoints/regionPoints/regionPoints.cpp
@@ -21,6 +22,8 @@ set(SourceFiles
     operation/postprocessOperation/postprocessOperation.cpp
     operation/PostprocessOperation/PostprocessOperationSum.cpp
     operation/PostprocessOperation/PostprocessOperationAverage.cpp
+    operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp
+    
     operation/includeMask/includeMask.cpp
     operation/includeMask/IncludeMasks.cpp
 
diff --git a/src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp b/src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp
index 508c572c..2ea4e81b 100644
--- a/src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp
+++ b/src/PostprocessData/fieldsDataBase/fieldsDataBase.cpp
@@ -21,17 +21,43 @@ Licence:
 #include <regex>
 
 #include "vocabs.hpp"
-#include "Time.hpp"
+#include "systemControl.hpp"
 #include "fieldsDataBase.hpp"
 #include "fieldFunctions.hpp"
-#include "dynamicPointStructure.hpp"
 
-bool pFlow::fieldsDataBase::checkForUpdate(const word &name, bool forceUpdate)
+
+namespace pFlow
+{
+
+bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
+
+}
+
+bool pFlow::fieldsDataBase::loadPointStructureToTime()
+{
+    return false;
+}
+
+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(name); found)
+    if(auto [iter, found]= captureTime_.findIf(compoundName); found)
     {
         shouldUpdate = iter->second < t || forceUpdate;
         iter->second = t;
@@ -39,7 +65,7 @@ bool pFlow::fieldsDataBase::checkForUpdate(const word &name, bool forceUpdate)
     else
     {
         shouldUpdate = true;
-        captureTime_.insertIf(name, t);
+        captureTime_.insertIf(compoundName, t);
     }
 
     return shouldUpdate;
@@ -70,6 +96,105 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetRealField(const word
         field.size());
 }
 
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetVolume(bool forceUpdate)
+{
+    const word fName = "volume";
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+        
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        auto vols = getShape().volume();
+        
+        FieldTypeHost<real> volField
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< volField.size(); i++)
+        {
+            volField[i] = vols[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(volField)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetDensity(bool forceUpdate)
+{
+    const word fName = "density";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        const auto index = updateFieldUint32("shapeIndex", true);
+        const auto dens = getShape().density();
+
+        FieldTypeHost<real> denFeild
+        (
+            fName, 
+            "value",
+            pointFieldSize()
+        );
+
+        for(uint32 i=0; i< denFeild.size(); i++)
+        {
+            denFeild[i] = dens[index[i]];
+        }
+
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            std::move(denFeild)  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
+pFlow::span<pFlow::real> pFlow::fieldsDataBase::createOrGetOne(bool forceUpdate)
+{
+    const word fName = "one";
+    
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
+
+    if(shouldUpdate)
+    {
+        allFields_.emplaceBackOrReplace<FieldTypeHost<real>>
+        (
+            fName,
+            FieldTypeHost<real> 
+            (
+                fName, 
+                "value",
+                pointFieldSize(),
+                1.0
+            )  
+        );
+    }
+
+    auto& field = allFields_.getObject<FieldTypeHost<real>>(fName);
+    return span<real>(
+        field.data(), 
+        field.size());
+}
+
 bool pFlow::fieldsDataBase::findFunction(
     const word &compoundFieldName,
     word &fieldName,
@@ -274,9 +399,9 @@ bool pFlow::fieldsDataBase::inputOutputType
     return false;
 }
 
-pFlow::fieldsDataBase::fieldsDataBase(const Time &time, bool inSimulation)
-: 
-    time_(time), 
+pFlow::fieldsDataBase::fieldsDataBase(systemControl& control, bool inSimulation)
+:
+    time_(control.time()),
     inSimulation_(inSimulation)
 {}
 
@@ -285,84 +410,110 @@ pFlow::timeValue pFlow::fieldsDataBase::currentTime() const
     return time_.currentTime();
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldTypeNameFunction
 (
-    const word& compoundName, 
-    word& fieldName,
+    const word& compoundName,
+    word& pointFieldName, 
     word& originalType, 
     word& typeAfterFunction,
     Functions& func
-)
+)const
 {
-    if( !findFunction(compoundName, fieldName, func))
+    if( !findFunction(compoundName, pointFieldName, func))
     {
-        fatalErrorInFunction;
+        fatalErrorInFunction
+            <<"Error in processing function for field: "
+            <<compoundName<<endl;;
         return false;
     }
 
-    if(reservedFieldNames_.contains(fieldName))
+    if( reservedFieldNames_.contains(pointFieldName))
     {
-        originalType = reservedFieldNames_.find(fieldName)->second;
+        // The name is in the reserved fields 
+        originalType = reservedFieldNames_.find(pointFieldName)->second;
     }
     else
     {
-        word fieldTypeName =  time_.lookupObjectTypeName(fieldName);
-        word space;
-        if(!pointFieldGetType(fieldTypeName, originalType, space))
+        // the name is in the Time object 
+        if(pointFieldNameExists(pointFieldName))
         {
-            fatalErrorInFunction<<"Cannot extract type name from "<< fieldTypeName<<endl;
+            originalType = getPointFieldType(pointFieldName);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "The field name: "<< pointFieldName
+                << " is not in the Time object.\n"
+                <<" Avaiable ones are: \n"<< time().objectNames()<<endl;
             return false;
         }
     }
-    
-    word outputType;
-    if(!inputOutputType(func, originalType, outputType))
+
+    if(!inputOutputType(func, originalType, typeAfterFunction))
     {
-        fatalErrorInFunction<<"Cannnot determine the input and output types for: "<< compoundName<<endl;
+        fatalErrorInFunction
+            <<"Cannnot determine the input and output types for: "
+            << compoundName<<endl;
         return false;
     }
-    typeAfterFunction = outputType;
+
     return true;
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldType
 (
-    const word &compoundName, 
-    word &originalType, 
-    word &typeAfterFunction
-)
+    const word& compoundName,
+    word& originalType, 
+    word& typeAfterFunction
+) const
 {
     Functions func;
     word fieldName;
-    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
 }
 
-bool pFlow::fieldsDataBase::getPointFieldType
+bool pFlow::fieldsDataBase::getFieldType
 (
     const word &compoundName, 
     word &typeAfterFunction
-)
+) const
 {
     Functions func;
     word originalType;
     word fieldName;
-    return getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func);
+    if( !getFieldTypeNameFunction(compoundName, fieldName, originalType, typeAfterFunction, func))
+    {
+        return false;
+    }    
+    return true;
 }
 
 pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updatePoints(bool forceUpdate)
 {
-    
-    bool shouldUpdate = checkForUpdate("position", forceUpdate);
+    const word fName = "position";
+    bool shouldUpdate = checkForUpdate(fName, forceUpdate);
         
     if(shouldUpdate)
     {
+        // load pointStructure 
+        if(!loadPointStructureToTime())
+        {
+            fatalErrorInFunction<< "Error in loading pointStructure to Time object."<<endl;
+            fatalExit;
+        }
         const auto& pstruct = pStruct(); 
-        allFields_.emplaceBackOrReplace<PointsTypeHost>(
-            "position", 
-            pstruct.activePointsHost());
+        allFields_.emplaceBackOrReplace<PointsTypeHost>
+        (
+            fName, 
+            pstruct.activePointsHost()
+        );
     }
    
-    auto& points = allFields_.getObject<PointsTypeHost>("position");
+    auto& points = allFields_.getObject<PointsTypeHost>(fName);
     
     return span<realx3>(
         points.data(), 
@@ -379,7 +530,12 @@ pFlow::span<pFlow::realx3> pFlow::fieldsDataBase::updateFieldRealx3
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func) )
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -410,7 +566,12 @@ pFlow::span<pFlow::realx4> pFlow::fieldsDataBase::updateFieldRealx4
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func))
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -441,7 +602,12 @@ pFlow::span<pFlow::real> pFlow::fieldsDataBase::updateFieldReal
     word originalType, typeAfterFunction, fieldName;
     Functions func;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldTypeNameFunction(
+        compoundName, 
+        fieldName, 
+        originalType, 
+        typeAfterFunction, 
+        func) )
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -609,11 +775,9 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
     bool forceUpdate
 )
 {
-
-    word originalType, typeAfterFunction, fieldName;
-    Functions func;
+    word typeAfterFunction, originalType;
     
-    if( !getPointFieldType(compoundName, fieldName, originalType, typeAfterFunction, func))
+    if( !getFieldType(compoundName, originalType, typeAfterFunction))
     {
         fatalErrorInFunction<< "Error in getting the type name of field: "<<
             compoundName<<", with type name: "<< originalType <<endl;
@@ -621,8 +785,7 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
         return span<real>(nullptr, 0);
     }
 
-
-    if(  typeAfterFunction== getTypeName<realx3>()  )
+    if(  typeAfterFunction == getTypeName<realx3>()  )
     {
         return updateFieldRealx3(compoundName, forceUpdate);
     }
@@ -642,19 +805,41 @@ pFlow::allPointFieldTypes pFlow::fieldsDataBase::updateFieldAll
     }
 }
 
-const pFlow::pointStructure &pFlow::fieldsDataBase::pStruct() const
+
+
+pFlow::uniquePtr<pFlow::fieldsDataBase> 
+    pFlow::fieldsDataBase::create(systemControl& control, bool inSimulation)
 {
-    if(inSimulation_)
+    word dbType;
+    if(inSimulation)
     {
-        return 
-            static_cast<const pointStructure&>(
-                time_.lookupObject<dynamicPointStructure>(pointStructureFile__)
-            );
+        dbType = "fieldsDataBase<simulation>";
     }
     else
     {
-        return time_.lookupObject<pointStructure>(pointStructureFile__);
+        dbType = "fieldsDataBase<postSimulation>";
     }
+    
+    if( boolvCtorSelector_.search(dbType) )
+    {
+        auto objPtr = 
+            boolvCtorSelector_[dbType](control, 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)
diff --git a/src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp b/src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp
index 13a71f31..ec8c577e 100644
--- a/src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp
+++ b/src/PostprocessData/fieldsDataBase/fieldsDataBase.hpp
@@ -21,52 +21,23 @@ Licence:
 #ifndef __fieldsDataBase_hpp__
 #define __fieldsDataBase_hpp__
 
-#include <variant>
-#include <concepts>
-
+#include "fieldsDataBaseDclr.hpp"
 
 #include "pointStructure.hpp"
 #include "pointFields.hpp"
 #include "anyList.hpp"
 #include "Map.hpp"
-#include "span.hpp"
+#include "shape.hpp"
+
 
 namespace pFlow
 {
 
+class systemControl;
 class Time;
 
-bool pointFieldGetType(const word& TYPENAME, word& fieldType, word& fieldSpace);
-
-template<typename T>
-concept ValidFieldType = 
-    std::same_as<T, real> || 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4> || 
-    std::same_as<T, uint32>;
-
-template<typename T> 
-concept VectorType = 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4>;
-
-template<typename T>
-concept ScalarType = 
-    std::same_as<T, real>;
-
-
-template<typename T>
-concept ValidRegionFieldType = 
-    std::same_as<T, real> || 
-    std::same_as<T, realx3> || 
-    std::same_as<T, realx4> ;
-
-
-using allPointFieldTypes = std::variant<span<real>, span<realx3>, span<realx4>>; 
-
 
 class fieldsDataBase
-
 {
 private:
 
@@ -99,9 +70,6 @@ private:
 
 
     // - Member variables
-    
-        /// const reference to the Time object
-        const Time&         time_;
 
         /// List to store all the point fields 
         anyList             allFields_;
@@ -109,114 +77,95 @@ private:
         /// 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"}
+            {"one", "real"},
+            {"volume", "real"},
+            {"density", "real"}
         };
-        
-    // - Private member functions
 
-        uint32 pointFieldSize()
-        {
-            auto s = updatePoints();
-            return s.size();
-        }
+        /// check if pointField name exists in Time or time folder 
+        virtual 
+        bool pointFieldNameExists(const word& name)const = 0;
         
-        /// @brief Checks if a field needs to be updated.
-        /// @param name The name of the field to check.
-        /// @param forceUpdate Forces an update if true. Defaults to `false`.
-        /// @return `true` if the field needs updating or `forceUpdate` is true, `false` otherwise.
-        bool checkForUpdate(const word& name, bool forceUpdate = false);
+        /// 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);
 
-        /**
-         * @brief Updates and retrieves a point field of a specified type from the database.
-         * 
-         * This is a template function that updates and retrieves a point field of type `T`
-         * from the database. It checks if the field needs to be updated based on the last
-         * capture time or if a forced update is requested. If an update is necessary, it
-         * retrieves the latest data for the field.
-         * 
-         * @tparam T The type of the point field to update and retrieve. Must be a ValidFieldType.
-         * @param name The name of the field to update.
-         * @param forceUpdate A boolean flag indicating whether to force an update of the field
-         *                    regardless of its current state. Defaults to `false`.
-         * 
-         * @return A span of `T` representing the updated field data.
-         * 
-         * @throws message If the field cannot be found in the database or if there is a type mismatch.
-         */
         template<ValidFieldType T>
         span<T> updateField(const word& name, bool forceUpdate = false);
+
+        template<ValidFieldType T>
+        span<T> updateReservedField(const word& name, bool forceUpdate = false);
         
-        /// @brief Creates a new real field or retrieves an existing one.
-        /// 
-        /// If a field with the given name already exists, it returns a span to that field.
-        /// If the field does not exist, it creates a new real field with the given name
-        /// and returns a span to the newly created field.
-        /// 
-        /// @param name The name of the field to create or retrieve.
-        /// @return span<real> of the field 
         span<real> createOrGetRealField(const word& name);
 
-        /**
-         * @brief Parses a compound field name to extract the base field name and function.
-         * 
-         * This function takes a compound field name, which may include a function applied
-         * to a base field (e.g., "mag(velocity)"), and parses it to extract the base field
-         * name (e.g., "velocity") and the function to be applied (e.g., `Functions::Magnitude`).
-         * 
-         * The function supports the following syntax for compound field names:
-         * - `fieldName` (no function applied)
-         * - `functionName(fieldName)`
-         * 
-         * Supported function names are defined in the `Functions` enum.
-         * 
-         * @param compoundFieldName The compound field name to parse.
-         * @param fieldName A reference to a `word` where the extracted base field name
-         *                  will be stored.
-         * @param func A reference to a `Functions` enum where the identified function
-         *             will be stored. If no function is applied, this will be set to
-         *             `Functions::None`.
-         * 
-         * @return `true` if the compound field name was successfully parsed and the base
-         *         field name and function were extracted, `false` otherwise.
-         * 
-         * @note The function modifies the `fieldName` and `func` parameters to return the
-         *       extracted information.
-         */
+        span<real> createOrGetVolume(bool forceUpdate=false);
+
+        span<real> createOrGetDensity(bool forceUpdate=false);
+
+        span<real> createOrGetOne(bool forceUpdate=false);
+
         static 
         bool findFunction(
             const word& compoundFieldName, 
             word& fieldName, 
             fieldsDataBase::Functions& func );
         
-        /**
-         * @brief Determines the input and output types for a given function.
-         * 
-         * This function takes a `Functions` enum value and an input type as a string
-         * and determines the corresponding output type based on the function being applied.
-         * 
-         * @param func The function for which to determine the input and output types.
-         * @param inputType The input type as a string.
-         * @param outputType A reference to a string where the determined output type will be stored.
-         * 
-         * @return `true` if the input and output types were successfully determined, `false` otherwise.
-         */
         static
         bool inputOutputType(
             fieldsDataBase::Functions func, 
             const word& inputType, 
             word& outputType);
 
+
+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(const Time& time, bool inSimulation);
+        fieldsDataBase(systemControl& control, bool inSimulation);
 
         /// no copy constructor
         fieldsDataBase(const fieldsDataBase&) = delete;
@@ -231,140 +180,70 @@ public:
         fieldsDataBase& operator=(fieldsDataBase&&) = delete;
         
         /// destructor
-        ~fieldsDataBase() = default;
+        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
+        /// const ref to object Time 
+        const Time& time()const 
         {
             return time_;
         }
-    
+        
+        Time& time()
+        {
+            return time_;
+        }
+
     // - Public Member Functions
 
-        /**
-         * @brief Retrieves the type of a point field based on its compound name.
-         * 
-         * This function attempts to extract the type of a point field from its compound name.
-         * The compound name may include additional information such as a function or operation
-         * applied to the field, ie. mag(velcoty). If the type is successfully determined, it  
-         * is stored in the provided `typeName` parameter.
-         * 
-         * @param compoundName The compound name of the field, which may include additional
-         *                     information about operations or functions applied to the field.
-         * @param originalType A reference to a `word` where the original type name is obtained.
-         *                     This will be set to the type of the field before any function is applied.
-         * @param typeAfterFunction A reference to a `word` where the type name after applying                   
-         *                          the function is obtained.
-         * @param func the applied function to the field.
-         *  
-         * @return `true` if the type was successfully determined and stored in `typeName`,
-         *         `false` otherwise.
-         */
-        bool getPointFieldType(
+        bool getFieldTypeNameFunction
+        (
             const word& compoundName,
-            word& fieldName, 
+            word& pointFieldName, 
             word& originalType, 
             word& typeAfterFunction,
-            Functions& func);   
-        
-        /// overload for the function getPointFieldType without `func` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& originalType, 
-            word& typeAfterFunction);
-        
-        /// overload for function getPointFieldType without `originalType` argument
-        bool getPointFieldType(
-            const word& compoundName, 
-            word& typeAfterFunction);
+            Functions& func
+        )const;
 
-        /**
-         * @brief Updates the points data and returns a span to the updated points.
-         * 
-         * This function ensures that the points data is up-to-date by checking if an update
-         * is necessary. If the data is outdated or if a forced update is requested, it retrieves
-         * the latest points data and stores it in the internal fields database. The function
-         * then returns a span to the updated points data for further use.
-         * 
-         * @param forceUpdate A boolean flag indicating whether to force an update of the points
-         *                    data regardless of its current state. Defaults to `false`.
-         * 
-         * @return A span of `realx3` representing the updated points data.
-         */
+        
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& originalType, 
+            word& typeAfterFunction
+        ) const;
+
+        bool getFieldType
+        (
+            const word& compoundName,
+            word& typeAfterFunction
+        ) const;
+
+        
+        /// update pointStructure if necessary 
         span<realx3> updatePoints(bool forceUpdate = false);
         
-        /**
-         * @brief Updates and retrieves a realx3 point field from the database.
-         * 
-         * This function retrieves or updates a realx3 field based on its compound name.
-         * The compound name cannot include any function operation.
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated realx3 field data
-         * 
-         * @throws message If the field type is not compatible with realx3 or if the field
-         *                 cannot be found in the database
-         */
+        /// update a field with realx3 type 
         span<realx3> updateFieldRealx3(
             const word& compoundName, 
             bool forceUpdate = false);
         
-        /**
-         * @brief Updates and retrieves a realx4 point field from the database.
-         * 
-         * This function retrieves or updates a realx4 field based on its compound name.
-         * The compound name cannot include any function operation.
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated realx3 field data
-         * 
-         * @throws message If the field type is not compatible with realx4 or if the field
-         *                 cannot be found in the database
-         */
         span<realx4> updateFieldRealx4(
             const word& compoundName, 
             bool forceUpdate = false);
         
-        /**
-         * @brief Updates and retrieves a real point field from the database.
-         * 
-         * This function retrieves or updates a real field based on its compound name.
-         * The compound name may include a function operation (e.g., abs, square, etc.).
-         * If the field needs to be updated or if forceUpdate is true, the method will
-         * fetch the latest data from the database and apply the specified function.
-         * 
-         * Supported functions include:
-         * - None: Retrieves the field as is.
-         * - abs: Computes the absolute value of the field.
-         * - square: Computes the square of the field.
-         * - cube: Computes the cube of the field.
-         * - sqrt: Computes the square root of the field.
-         * - mag, magSquare, magCube, magSqrt: Compute magnitude-related operations for vector fields.
-         * - component(x, y, z): Extracts a specific component from a vector field.
-         * 
-         * @param compoundName The name of the field, possibly including a function operation.
-         * @param forceUpdate If true, forces an update of the field regardless of its current state.
-         *                    Defaults to false.
-         * 
-         * @return A span containing the updated real field data.
-         * 
-         * @throws message If the field type is not compatible with real or if the field
-         *                 cannot be found in the database.
-         */
         span<real> updateFieldReal(
             const word& compoundName, 
             bool forceUpdate = false);
@@ -378,11 +257,20 @@ public:
         allPointFieldTypes updateFieldAll(
             const word& compoundName, 
             bool forceUpdate = false);
+        
+        virtual 
+        const pointStructure& pStruct()const = 0;
 
-        const pointStructure& pStruct()const;
+        virtual 
+        void resetTimeFolder() = 0;
+    
+    static
+    uniquePtr<fieldsDataBase> create(
+        systemControl& control, 
+        bool inSimulation);
 };
 
-} // nampespace pFlow
+} // namespace pFlow
 
 #include "fieldsDataBaseTemplates.cpp"
 
diff --git a/src/PostprocessData/fieldsDataBase/fieldsDataBaseDclr.hpp b/src/PostprocessData/fieldsDataBase/fieldsDataBaseDclr.hpp
new file mode 100644
index 00000000..283b5e7b
--- /dev/null
+++ b/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__
\ No newline at end of file
diff --git a/src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp b/src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp
index 9466a8d6..a36eb95c 100644
--- a/src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp
+++ b/src/PostprocessData/fieldsDataBase/fieldsDataBaseTemplates.cpp
@@ -32,40 +32,25 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
         
     if(shouldUpdate)
     {
-        if(name == "one")
-        {   
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>
-            (
-                "one",
-                FieldTypeHost<T>
-                (
-                    "one", 
-                    "value",
-                    pointFieldSize(),
-                    T{1}
-                )
-            );
-        }
-        else if( name == "position")
+        if(reservedFieldNames_.contains(name))
         {
-            if constexpr( std::same_as<T, realx3>)
+            return updateReservedField<T>(name, true);
+        }   
+        else
+        {
+            if( loadPointFieldToTime(name) )
             {
-                return updatePoints(true);            
+                const auto& pField = time_.lookupObject<pointField_D<T>>(name);
+                allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
+                name, 
+                pField.activeValuesHost());
             }
             else
             {
-                fatalErrorInFunction<< "Error in getting the type name of field: "<<
-                    name<<", with type name: "<< getTypeName<T>() <<endl;
+                fatalErrorInFunction<<"Error in loading the pointField "<<name<<endl;
                 fatalExit;
-                return span<T>(nullptr, 0);
             }
-        }
-        else
-        {
-            const auto& pField = time_.lookupObject<pointField_D<T>>(name);
-            allFields_.emplaceBackOrReplace<FieldTypeHost<T>>(
-            name, 
-            pField.activeValuesHost());
+            
         }
     }
 
@@ -77,4 +62,80 @@ pFlow::span<T> pFlow::fieldsDataBase::updateField(const word& name, bool forceUp
     
 }
 
+
+template<pFlow::ValidFieldType T>
+inline
+pFlow::span<T> pFlow::fieldsDataBase::updateReservedField
+(
+    const word& name, 
+    bool forceUpdate
+)
+{
+    if(name == "one")
+    {   
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetOne(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field one."<<endl;
+            fatalExit;
+        }
+    }
+    else if(name == "volume")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetVolume(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field volume."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "density")
+    {
+        if constexpr( std::same_as<T,real>)
+        {
+            return createOrGetDensity(forceUpdate);
+        }
+        else
+        {
+            fatalErrorInFunction
+                << "This type: "
+                << getTypeName<T>()
+                <<" is not supported for field density."<<endl;
+            fatalExit;
+        }
+    }
+    else if( name == "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__
\ No newline at end of file
diff --git a/src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.cpp b/src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.cpp
new file mode 100644
index 00000000..312f8e51
--- /dev/null
+++ b/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()
+{
+}
diff --git a/src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.hpp b/src/PostprocessData/fieldsDataBase/simulationFieldsDataBase.hpp
new file mode 100644
index 00000000..d7ea7f50
--- /dev/null
+++ b/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__
\ No newline at end of file
diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp
index 16b029f1..5c0f164e 100644
--- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp
+++ b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp
@@ -128,7 +128,7 @@ Licence:
 #include <variant>
 #include <vector>
 
-#include "postprocessOperation.hpp"
+#include "PostprocessOperationAverage.hpp"
 #include "regionField.hpp"
 #include "includeMask.hpp"
 
@@ -138,7 +138,7 @@ namespace pFlow
 
 class PostprocessOperationAvMassVelocity
 :
-    public postprocessOperation
+    public PostprocessOperationAverage
 {
    
 public:
diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp
index a0de2384..3ffbd5b8 100644
--- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp
+++ b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp
@@ -1,7 +1,7 @@
 #include "PostprocessOperationAverage.hpp"
 #include "dictionary.hpp"
 #include "fieldsDataBase.hpp"
-#include "fieldFunctions.hpp"
+#include "operationFunctions.hpp"
 
 /// Constructs average processor and initializes result field based on input field type
 pFlow::PostprocessOperationAverage::PostprocessOperationAverage
@@ -80,7 +80,8 @@ pFlow::PostprocessOperationAverage::PostprocessOperationAverage
 /// Performs weighted average of field values within each region
 bool pFlow::PostprocessOperationAverage::execute
 (
-    const std::vector<span<real>>& weights
+    const std::vector<span<real>>& weights,
+    const regionField<real>& volFactor
 )
 {
     auto allField = database().updateFieldAll(fieldName());
@@ -99,7 +100,7 @@ bool pFlow::PostprocessOperationAverage::execute
                 return executeAverageOperation(
                 procName,
                 field,
-                regP,
+                volFactor,
                 dbVol,
                 weights,
                 phi,
@@ -124,6 +125,7 @@ bool pFlow::PostprocessOperationAverage::execute
                         procName,
                         field,
                         std::get<regionField<T>>(processedRegField),
+                        volFactor,
                         dbVol,
                         weights,
                         mask);
diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp
index 93181c1e..6620d588 100644
--- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp
+++ b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.hpp
@@ -193,7 +193,9 @@ public:
     /// @brief Execute average operation on field values
     /// @param weights Weight factors for particles
     /// @return True if successful
-    bool execute(const std::vector<span<real>>& weights) override;
+    bool execute(
+        const std::vector<span<real>>& weights,
+        const regionField<real>& volFactor) override;
 
 }; 
 
diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp
index 922bc1d4..c10315d7 100644
--- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp
+++ b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp
@@ -1,7 +1,7 @@
 #include "PostprocessOperationSum.hpp"
 #include "dictionary.hpp"
 #include "fieldsDataBase.hpp"
-#include "fieldFunctions.hpp"
+#include "operationFunctions.hpp"
 
 /// Constructs sum processor and initializes result field based on input field type
 pFlow::PostprocessOperationSum::PostprocessOperationSum
@@ -41,7 +41,8 @@ pFlow::PostprocessOperationSum::PostprocessOperationSum
 /// Performs weighted sum of field values within each region
 bool pFlow::PostprocessOperationSum::execute
 (
-    const std::vector<span<real>>& weights
+    const std::vector<span<real>>& weights,
+    const regionField<real>& volFactor
 )
 {
     auto allField = database().updateFieldAll(fieldName());
@@ -60,7 +61,7 @@ bool pFlow::PostprocessOperationSum::execute
                 return executeSumOperation(
                 procName,
                 field,
-                regP,
+                volFactor,
                 dbVol,
                 weights,
                 phi,
diff --git a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp
index c8e910ff..b811e89b 100644
--- a/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp
+++ b/src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp
@@ -175,7 +175,9 @@ public:
     /// @brief Execute sum operation on field values
     /// @param weights Weight factors for particles
     /// @return True if successful
-    bool execute(const std::vector<span<real>>& weights) override;
+    bool execute(
+        const std::vector<span<real>>& weights,
+        const regionField<real>& volFactor) override;
 
 }; 
 
diff --git a/src/PostprocessData/operation/PostprocessOperation/fieldFunctions.hpp b/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp
similarity index 88%
rename from src/PostprocessData/operation/PostprocessOperation/fieldFunctions.hpp
rename to src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp
index 3e278858..8dbe643b 100644
--- a/src/PostprocessData/operation/PostprocessOperation/fieldFunctions.hpp
+++ b/src/PostprocessData/operation/PostprocessOperation/operationFunctions.hpp
@@ -18,8 +18,8 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-#ifndef __fieldFunctions_hpp__
-#define __fieldFunctions_hpp__
+#ifndef __operationFunctions_hpp__
+#define __operationFunctions_hpp__
 
 #include <vector>
 
@@ -36,13 +36,14 @@ regionField<T> executeSumOperation
 (
     const word&                     regFieldName,
     const span<T>&                  field, 
-    const regionPoints&             regPoints,
+    const regionField<real>&        volFactor,
     const bool                      devideByVol,
     const std::vector<span<real>>&  weights,
     const span<real>&               phi,
     const includeMask::Mask&        mask
 )
 {
+    const auto& regPoints = volFactor.regPoints(); 
     regionField<T> processedField(regFieldName, regPoints, T{});
     auto vols = regPoints.volumes();
 
@@ -63,7 +64,7 @@ regionField<T> executeSumOperation
         }
         if(devideByVol)
         {
-            processedField[reg] = sum/vols[reg];
+            processedField[reg] = sum/(volFactor[reg] * vols[reg]);
         }
         else
         {
@@ -80,13 +81,15 @@ regionField<T> executeAverageOperation
 (
     const word&                    regFieldName,
     const span<T>&                 field, 
-    const regionPoints&            regPoints,
-    const bool                     devideByVol, 
+    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();
 
@@ -113,7 +116,7 @@ regionField<T> executeAverageOperation
         
         if(devideByVol)
         {
-            processedField[reg] = sumNum / max(sumDen, smallValue) / vols[reg];
+            processedField[reg] = sumNum / max(sumDen, smallValue) / (volFactor[reg] * vols[reg]);
         }
         else
         {
@@ -131,7 +134,8 @@ regionField<T> executeFluctuation2Operation
 (
     const word&                    regFieldName,
     const span<T>&                 field,
-    const regionField<T>&          fieldAvg, 
+    const regionField<T>&          fieldAvg,
+    const regionField<real>&       volFactor, 
     const bool                     devideByVol, 
     const std::vector<span<real>>& weights,
     const includeMask::Mask&       mask
@@ -145,7 +149,7 @@ regionField<T> executeFluctuation2Operation
     {
         auto partIndices = regPoints.indices(reg);
         auto w = weights[reg];
-        auto vol = vols[reg];
+        auto vol = volFactor[reg] * vols[reg];
         T avField{};
         if(devideByVol)
         {
@@ -188,4 +192,4 @@ regionField<T> executeFluctuation2Operation
 
 } // namespace pFlow
 
-#endif //__fieldFunctions_hpp__
+#endif //__operationFunctions_hpp__
diff --git a/src/PostprocessData/operation/includeMask/IncludeMask.hpp b/src/PostprocessData/operation/includeMask/IncludeMask.hpp
index 6a4a24f7..34ae4dfb 100644
--- a/src/PostprocessData/operation/includeMask/IncludeMask.hpp
+++ b/src/PostprocessData/operation/includeMask/IncludeMask.hpp
@@ -154,11 +154,34 @@ public:
             ).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(
+	add_vCtor
+    (
 		includeMask,
 		IncludeMask,
-		dictionary);
+		dictionary
+    );
+
+    add_vCtor
+    (
+		includeMask,
+		IncludeMask,
+		word
+    );
 
     /// Returns the mask for filtering elements (updates the mask if necessary)
     Mask getMask() override
@@ -203,20 +226,40 @@ public:
 
 	TypeInfoTemplate12("IncludeMask", T, allOp<int8>);
 
-	IncludeMask(
-		const dictionary& opDict, 
+    IncludeMask(
+        const dictionary& opDict, 
         fieldsDataBase& feildsDB)
-	:
-		includeMask(opDict, feildsDB)
-	{
+    :
+        includeMask(opDict, feildsDB)
+    {
         span<realx3>  s = database().updatePoints();        
         mask_.resize(s.size(), true);
-	}
+    }
 
-	add_vCtor(
-		includeMask,
-		IncludeMask,
-		dictionary);
+    IncludeMask(
+        const word& type,
+        const dictionary& opDict, 
+        fieldsDataBase& feildsDB)
+    :
+        includeMask(type, opDict, feildsDB)
+    {
+        span<realx3>  s = database().updatePoints();        
+        mask_.resize(s.size(), true);
+    }
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        dictionary
+    );
+
+    add_vCtor
+    (
+        includeMask,
+        IncludeMask,
+        word
+    );
 
     Mask getMask()override
     {
diff --git a/src/PostprocessData/operation/includeMask/includeMask.cpp b/src/PostprocessData/operation/includeMask/includeMask.cpp
index 95c85fb8..209391be 100644
--- a/src/PostprocessData/operation/includeMask/includeMask.cpp
+++ b/src/PostprocessData/operation/includeMask/includeMask.cpp
@@ -57,7 +57,7 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create
 		auto& maskDict = opDict.subDict(mask+"Info");
 		word maskField = maskDict.getVal<word>("field");
 		
-		if( !fieldsDB.getPointFieldType(maskField, fieldType) )
+		if( !fieldsDB.getFieldType(maskField, fieldType) )
 		{
 			fatalErrorInFunction<<"Error in retriving the type of field"
 				<< maskField <<" from dictionary "
@@ -111,7 +111,7 @@ pFlow::uniquePtr<pFlow::includeMask>
 		auto& maskDict = opDict.subDict(type+"Info");
 		word maskField = maskDict.getVal<word>("field");
 		
-		if( !fieldsDB.getPointFieldType(maskField, fieldType) )
+		if( !fieldsDB.getFieldType(maskField, fieldType) )
 		{
 			fatalErrorInFunction<<"Error in retriving the type of field"
 				<< maskField <<" from dictionary "
@@ -143,7 +143,7 @@ pFlow::uniquePtr<pFlow::includeMask>
 			method << " dose not exist. \n"
 			<<"Avaiable ones are: \n\n"
 			,
-			dictionaryvCtorSelector_
+			wordvCtorSelector_
 		);
 		fatalExit;
 		return nullptr;
diff --git a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp
index 53473a09..e752f359 100644
--- a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp
+++ b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp
@@ -117,7 +117,7 @@ pFlow::postprocessOperation::postprocessOperation
     ),
     divideByVolume_
     (
-        opDict.getValOrSet<Logical>("dividedByVolume", Logical(false))
+        opDict.getValOrSet<Logical>("divideByVolume", Logical(false))
     ),
     regionPoints_
     (
@@ -141,7 +141,7 @@ pFlow::postprocessOperation::postprocessOperation
     )
 {
 
-    if(!fieldsDB.getPointFieldType(fieldName_, fieldType_))
+    if(!fieldsDB.getFieldType(fieldName_, fieldType_))
     {
         fatalErrorInFunction;
         fatalExit;
@@ -159,7 +159,7 @@ bool pFlow::postprocessOperation::write(const fileSystem &parDir) const
     if(!osPtr_)
     {
         fileSystem path = parDir+(
-            processedFieldName() + ".Start_" + ti.prevTimeName());
+            processedFieldName() + ".Start_" + ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
         
         regPoints().write(osPtr_());
diff --git a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp
index 86e776d3..285a7b7e 100644
--- a/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp
+++ b/src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp
@@ -246,7 +246,10 @@ public:
 
     /// execute the operation
     /// @param weights Vector of weights for the operation.
-    virtual bool execute(const std::vector<span<real>>& weights) = 0;
+    /// @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.
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp
index 759395f5..10590f7c 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp
@@ -36,6 +36,12 @@ pFlow::PostprocessComponent<RegionType,ProcessMethodType>::PostprocessComponent
     (
         regionPointsPtr_().size()
     ),
+    volumeFactor_
+    (
+        "volumeFactor",
+        regionPointsPtr_(),
+        1.0
+    ),
     operationDicts_(readDictList("operations", dict))
 {
 
@@ -106,7 +112,7 @@ bool pFlow::PostprocessComponent<RegionType, ProcessMethodType>::execute
     
     for(auto& op:operatios_)
     {
-        if( !op->execute(weights) )
+        if( !op->execute(weights, volumeFactor_) )
         {
             fatalErrorInFunction
                 <<"error occured in executing operatoin defined in dict "
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp
index 51c46490..996467a1 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.hpp
@@ -51,6 +51,8 @@ private:
     /// Method for processing the selected particles data 
     std::vector<ProcessMethodType>  regionsProcessMethod_;
 
+    regionField<real>               volumeFactor_;
+
     bool    executed_{false};
 
     dictionaryList                  operationDicts_;
@@ -62,6 +64,16 @@ protected:
         return regionsProcessMethod_;
     }
 
+    regionField<real>& volumeFactor()
+    {
+        return volumeFactor_;
+    }
+    
+    const regionField<real>& volumeFactor()const
+    {
+        return volumeFactor_;
+    }
+
 public:
 
     // type info 
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp
index d801c9d4..7b9a37c3 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentArithmetic.hpp
@@ -51,10 +51,12 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
             regs[i] = arithmetic(); // Changed from uniformDistribution() to arithmetic()
+            volFactor[i] = 1.0;
         }
     }
 
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp
index d846ddcc..0fb00163 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentGaussian.hpp
@@ -23,6 +23,7 @@ Licence:
 
 #include "PostprocessComponent.hpp"
 #include "GaussianDistribution.hpp"
+#include "numericConstants.hpp"
 
 namespace pFlow
 {
@@ -51,15 +52,20 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
+
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
-            regs[i] = GaussianDistribution(c[i], pow(d[i],2));
+            auto r = d[i]/2;
+            regs[i] = GaussianDistribution(c[i], pow(r/3.0,2));
+            volFactor[i] = 0.677683 / (4.0/3.0*Pi*r);
         }
     }
 
     // add the virtual constructor
-    add_vCtor(
+    add_vCtor
+    (
         postprocessComponent,
         PostprocessComponentGaussian,
         dictionary
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp
index 29f0fda9..5a181474 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponentUniform.hpp
@@ -51,10 +51,12 @@ public:
         auto d = this->regPoints().eqDiameters();
         auto c = this->regPoints().centers();
         auto& regs = this->regionProecessMethod();
+        auto& volFactor = this->volumeFactor();
         const uint32 n = d.size();
         for(uint32 i=0; i<n; i++)
         {
             regs[i] = uniformDistribution();
+            volFactor[i] = 1.0;
         }
     }
 
diff --git a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp
index 71bcb054..2a6245a5 100644
--- a/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp
+++ b/src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponents.cpp
@@ -24,10 +24,19 @@ Licence:
 
 // region types 
 #include "sphereRegionPoints.hpp"
-
+#include "lineRegionPoints.hpp"
+#include "multipleSpheresRegionPoints.hpp"
 
 
 template class pFlow::PostprocessComponentGaussian<pFlow::sphereRegionPoints>;
 template class pFlow::PostprocessComponentUniform<pFlow::sphereRegionPoints>;
 template class pFlow::PostprocessComponentArithmetic<pFlow::sphereRegionPoints>;
 
+template class pFlow::PostprocessComponentGaussian<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::multipleSpheresRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::multipleSpheresRegionPoints>;
+
+template class pFlow::PostprocessComponentGaussian<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentUniform<pFlow::lineRegionPoints>;
+template class pFlow::PostprocessComponentArithmetic<pFlow::lineRegionPoints>;
+
diff --git a/src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp b/src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp
index f388c939..7588e602 100644
--- a/src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp
+++ b/src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp
@@ -135,7 +135,7 @@ bool pFlow::particleProbePostprocessComponent::write(const fileSystem& parDir)co
     if( !osPtr_)
     {
         // file is not open yet
-        fileSystem path = parDir + (name_+".Start_"+ti.prevTimeName());
+        fileSystem path = parDir + (name_+".Start_"+ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
         regionPointsPtr_().write(osPtr_());
     }
diff --git a/src/PostprocessData/postprocessData/postprocessData.cpp b/src/PostprocessData/postprocessData/postprocessData.cpp
index 15a7709b..dc1dc479 100644
--- a/src/PostprocessData/postprocessData/postprocessData.cpp
+++ b/src/PostprocessData/postprocessData/postprocessData.cpp
@@ -30,7 +30,14 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
 :
     auxFunctions(control),
     time_(control.time()),
-    fieldsDataBase_(control.time(), true),
+    fieldsDataBasePtr_
+    ( 
+        fieldsDataBase::create
+        (
+            const_cast<systemControl&>(control),
+            true
+        )
+    ),
     dict_
     (
         objectFile
@@ -73,6 +80,12 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
             control.time().saveInterval(),
             "execution");
     }
+
+    shapeType_ = dict_.getValOrSet<word>
+    (
+        "shapeType", 
+        word("sphere")
+    );
     
     componentsDictsPtr_ = makeUnique<dictionaryList>(readDictList("components", dict_));
 
@@ -80,7 +93,7 @@ pFlow::postprocessData::postprocessData(const systemControl &control)
     {
         postprocesses_.push_back( postprocessComponent::create(
             compDict, 
-            fieldsDataBase_, 
+            fieldsDataBasePtr_(), 
             defaultTimeControlPtr_() ));
     }
 
diff --git a/src/PostprocessData/postprocessData/postprocessData.hpp b/src/PostprocessData/postprocessData/postprocessData.hpp
index 4e5a95f2..e2223116 100644
--- a/src/PostprocessData/postprocessData/postprocessData.hpp
+++ b/src/PostprocessData/postprocessData/postprocessData.hpp
@@ -57,11 +57,14 @@ class postprocessData
     const Time&           time_;
 
     /// Database for all the points fields on the host 
-    fieldsDataBase        fieldsDataBase_;
+    uniquePtr<fieldsDataBase>        fieldsDataBasePtr_;
 
     /// file dictionary that is constructed from the file (postProcessDataDict)
     fileDictionary        dict_;
 
+    /// 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;
 
diff --git a/src/PostprocessData/postprocessData/postprocessGlobals.hpp b/src/PostprocessData/postprocessData/postprocessGlobals.hpp
index 99d04038..335bc267 100644
--- a/src/PostprocessData/postprocessData/postprocessGlobals.hpp
+++ b/src/PostprocessData/postprocessData/postprocessGlobals.hpp
@@ -27,7 +27,7 @@ namespace pFlow::postProcessGlobals
 {
 
 static fileSystem defaultDir__;
-inline const word defaultRelDir__ = "VTK/postprocessData";
+inline const word defaultRelDir__ = "postprocessData";
 
 }
 
diff --git a/src/PostprocessData/processMethod/GaussianDistribution.hpp b/src/PostprocessData/processMethod/GaussianDistribution.hpp
index 53889b76..6248bfa0 100644
--- a/src/PostprocessData/processMethod/GaussianDistribution.hpp
+++ b/src/PostprocessData/processMethod/GaussianDistribution.hpp
@@ -26,7 +26,7 @@ Licence:
 #include "typeInfo.hpp"
 #include "types.hpp"
 #include "span.hpp"
-
+#include "numericConstants.hpp"
 
 namespace pFlow
 {
@@ -74,15 +74,15 @@ public:
         for(uint32 i=0; i<indices.size(); i++)
         {
             auto x = points[indices[i]]-center;
-            auto f = exp(- dot(x,x)/(2*variance_));
+            auto f = exp(- dot(x,x)/(2*variance_))/sqrt(2.0*Pi*variance_);
             weight_[i] = f;
             sum += f;
         }
 
-        for(auto& w: weight_)
+        /*for(auto& w: weight_)
         {
             w /= sum;
-        }   
+        } */  
         return true;
     }
 
diff --git a/src/PostprocessData/region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp b/src/PostprocessData/region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
index e16ef774..b5ffb730 100644
--- a/src/PostprocessData/region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
+++ b/src/PostprocessData/region/regionPoints/centerPointsRegionPoints/centerPointsRegionPoints.cpp
@@ -23,11 +23,13 @@ bool pFlow::centerPointsRegionPoints::selectIds()
         }
     }
     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_);
+            probDict_.subDict(selector+"Info"));
         auto selectedPoints = selectorPtr->selectedPoints();
         ids_.resize(selectedPoints.size());
         ids_.assign(selectedPoints.begin(), selectedPoints.end());
diff --git a/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.cpp b/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
index 5e41ddf3..979a51a6 100644
--- a/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
+++ b/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.cpp
@@ -38,6 +38,7 @@ pFlow::lineRegionPoints::lineRegionPoints
     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)
@@ -64,6 +65,21 @@ pFlow::span<const pFlow::uint32> pFlow::lineRegionPoints::indices(uint32 elem) c
         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();
diff --git a/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.hpp b/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.hpp
index bdffd094..58989aab 100644
--- a/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.hpp
+++ b/src/PostprocessData/region/regionPoints/lineRegionPoints/lineRegionPoints.hpp
@@ -126,6 +126,8 @@ public:
 
     /// 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;
diff --git a/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp b/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
index 8665f72e..e983c602 100644
--- a/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
+++ b/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.cpp
@@ -59,6 +59,20 @@ pFlow::span<const pFlow::uint32> pFlow::multipleSpheresRegionPoints::indices(uin
     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();
diff --git a/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.hpp b/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.hpp
index f250d216..e01da658 100644
--- a/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.hpp
+++ b/src/PostprocessData/region/regionPoints/multipleSpheresRegionPoints/multipleSpheresRegionPoints.hpp
@@ -138,6 +138,11 @@ public:
     /// @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;
diff --git a/src/PostprocessData/sampleDictionary/postprocessDataDict b/src/PostprocessData/sampleDictionary/postprocessDataDict
index 3cf582a1..92e7d8a5 100755
--- a/src/PostprocessData/sampleDictionary/postprocessDataDict
+++ b/src/PostprocessData/sampleDictionary/postprocessDataDict
@@ -50,7 +50,7 @@ components
               {
                   function	  average;
                   field           velocity;
-                  dividedByVolume 	  no;  //default 
+                  divideByVolume 	  no;  //default 
                   threshold 	  3;   //default is 1;
                   includeMask     all;
               }
@@ -61,7 +61,7 @@ components
                   function        average;
                   field           one;
                   phi             one; // default
-                  dividedByVolume no;
+                  divideByVolume  no;
                   includeMask     lessThan;
                   
                   // diameter of par1 is 0.003, so these settings 
@@ -79,7 +79,7 @@ components
                   function        sum;
                   field           one;
                   phi             one; // default
-                  dividedByVolume yes;
+                  divideByVolume  yes;
                   
               }
           );
@@ -114,14 +114,14 @@ components
               {
                   function	   sum;
                   field            one;
-                  dividedByVolume  yes;  //default is no 
+                  divideByVolume  yes;  //default is no 
               }
               
               volumeDensity
               {
                 function         sum;
                 field            cube(diameter); // d^3, although it differs by pi/6
-                dividedByVolume  yes;  //default is no 
+                divideByVolume   yes;  //default is no 
               }
          );
     }
diff --git a/src/phasicFlow/repository/repository/repository.cpp b/src/phasicFlow/repository/repository/repository.cpp
index 648b912f..2815956f 100644
--- a/src/phasicFlow/repository/repository/repository.cpp
+++ b/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 )
diff --git a/src/phasicFlow/repository/repository/repository.hpp b/src/phasicFlow/repository/repository/repository.hpp
index 5e712c16..43984dc0 100644
--- a/src/phasicFlow/repository/repository/repository.hpp
+++ b/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);
diff --git a/src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp b/src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp
index 3a4a28f0..2e1591de 100644
--- a/src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp
+++ b/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()));