bug fix for pFlowToVTK for issue #72 in Plus branch

fixed the error when insert particles from file with integer number

benchmarks/helicalMixer/helicalMixer_1m/settings/settingsDict

@@ -12,7 +12,7 @@ dt                  0.00001;      // time step for integration (s)
 
 startTime           0;            // start time for simulation
 
-endTime             2;          // end time for simulation
+endTime             7.5;          // end time for simulation
 
 saveInterval        0.05;         // time interval for saving the simulation
 

benchmarks/helicalMixer/helicalMixer_2m/settings/settingsDict

@@ -12,7 +12,7 @@ dt                  0.00001;      // time step for integration (s)
 
 startTime           0;            // start time for simulation
 
-endTime             2;          // end time for simulation
+endTime             7.5;          // end time for simulation
 
 saveInterval        0.05;         // time interval for saving the simulation
 

benchmarks/helicalMixer/helicalMixer_4m/settings/settingsDict

@@ -12,7 +12,7 @@ dt                  0.00001;      // time step for integration (s)
 
 startTime           0;            // start time for simulation
 
-endTime             2;          // end time for simulation
+endTime             7.5;          // end time for simulation
 
 saveInterval        0.05;         // time interval for saving the simulation
 

benchmarks/helicalMixer/helicalMixer_500k/settings/settingsDict

@@ -12,7 +12,7 @@ dt                  0.00001;      // time step for integration (s)
 
 startTime           0;            // start time for simulation
 
-endTime             2;          // end time for simulation
+endTime             7.5;          // end time for simulation
 
 saveInterval        0.05;         // time interval for saving the simulation
 

benchmarks/helicalMixer/readme.md

@@ -1 +1,101 @@
 # Helical Mixer Benchmark (phasicFlow v-1.0)
+
+## Overview
+
+This benchmark compares the performance of phasicFlow with a well-stablished commercial DEM software for simulating a helical mixer with varying particle counts (250k to 4M particles). The benchmark measures both computational efficiency and memory usage across different hardware configurations.
+
+**Summary of Results:**
+
+- phasicFlow achieves similar performance to the commercial DEM software on the same hardware.
+- phasicFlow shows a 30% performance improvement when using the NVIDIA RTX A4000 compared to the RTX 4050Ti.
+- Memory usage is approximately 50% lower in phasicFlow compared to the commercial software, with phasicFlow using about 0.7 GB of memory per million particles, while the commercial software uses about 1.5 GB per million particles.
+
+## Simulation Setup
+
+<div align="center">
+    <img src="./images/commericalDEMsnapshot.png"/>
+    <div align="center">  
+        <p>Figure 1. Commercial DEM simulation snapshot</p>
+    </div>
+</div>
+
+<div align="center">
+    <img src="./images/phasicFlow_snapshot.png"/>
+    <div align="center">  
+        <p>Figure 2. phasicFlow simulation snapshot and visualized using Paraview</p>
+    </div>
+</div>
+
+### Hardware Specifications
+
+<div align="center">
+    Table 1. Hardware specifications used for benchmarking.
+</div>
+
+|  System     |           CPU            |             GPU              | Operating System |
+| :---------: | :----------------------: | :--------------------------: | :--------------: |
+|   Laptop    | Intel i9-13900HX 2.2 GHz | NVIDIA GeForce RTX 4050Ti 6G | Windows 11 24H2  |
+| Workstation | Intel Xeon 4210 2.2 GHz  |     NVIDIA RTX A4000 16G     |   Ubuntu 22.04   |
+
+### Simulation Parameters
+
+<div align="center">
+    Table 2. Parameters for helical mixer simulations.
+</div>
+
+| Case     | Particle Diameter | Particle Count | 
+| :-------: | :---------------: | :--------------: | 
+| 250k     | 6 mm              | 250,000        | 
+| 500k     | 5 mm              | 500,000        | 
+| 1M       | 4 mm              | 1,000,000      | 
+| 2M       | 3 mm              | 2,000,000      | 
+| 4M       | 2 mm              | 4,000,000      |
+
+The time step for all simulations was set to 1.0e-5 seconds and the simulation ran for 7.5 seconds.
+
+## Performance Comparison
+
+### Execution Time
+
+<div align="center">
+    Table 3. Total calculation time (minutes) for different configurations.
+</div>
+
+|     Software      | 250k   | 500k   | 1M     | 2M     | 4M     | 
+| :---------------: | :----: | :-----: | :-----: | :-----: | :-----: |
+| phasicFlow-4050Ti | 110 min | 215 min | 413 min | - | -       | 
+| Commercial DEM-4050Ti | 111 min | 210 min | 415 min | - | -       | 
+| phasicFlow-A4000  | 82 min | 150 min  | 300 min | 613 min | 1236 min |
+
+The execution time scales linearly with particle count. phasicFlow demonstrates approximately:
+
+- the computing speed is basically the same as well-established commercial DEM software on the same hardware
+- 30% performance improvement when using the NVIDIA RTX A4000 compared to the RTX 4050Ti
+
+<div align="center">
+    <img src="./images/performance.png"/>
+    <p>Figure 3. Calculation time comparison between phasicFlow and the well-established commercial DEM software.</p>
+</div>
+
+### Memory Usage
+
+<div align="center">
+    Table 4. Memory consumption for different configurations.
+</div>
+
+|     Software      | 250k   | 500k   | 1M      | 2M      | 4M      | 
+| :---------------: | :-----: | :-----: | :-----: | :-----: | :-----: | 
+| phasicFlow-4050Ti | 260 MB  | 404 MB  | 710 MB  | - | -       | 
+| Commercial DEM-4050Ti | 460 MB | 920 MB | 1574 MB | - | -       | 
+| phasicFlow-A4000  | 352 MB  | 496 MB  | 802 MB  | 1376 MB | 2310 MB |
+
+Memory efficiency comparison:
+
+- phasicFlow uses approximately 0.7 GB of memory per million particles
+- Commercial DEM software uses approximately 1.5 GB of memory per million particles
+- phasicFlow shows ~50% lower memory consumption compared to the commercial alternative
+- The memory usage scales linearly with particle count in both software packages. But due to memory limitations on GPUs, it is possible to run larger simulation on GPUs with phasicFlow.
+
+## Run Your Own Benchmarks
+
+The simulation case setup files are available in this folder for users interested in performing similar benchmarks on their own hardware. These files can be used to reproduce the tests and compare performance across different systems.

benchmarks/readme.md

@@ -1,7 +1,9 @@
 
 # Benchmarks
 
-Benchmakrs has been done on two different simulations: a simulation with simple geometry (rotating drum) and a simulation with complex geometry (helical mixer).
+Benchmakrs has been done on two different simulations: simulation with simple geometry (rotating drum) and a simulation with complex geometry (helical mixer). These benchmarks are used to show how PhasicFlow performs in different scenarios.
 
-- [rotating drum](./rotatingDrum/readme.md)
-- [helical mixer](./helicalMixer/readme.md)
+- [rotating drum](./rotatingDrum/)
+- [helical mixer](./helicalMixer/)
+
+**Note:** If you have performed benchmarks with PhasicFlow using other hardware or software other than PhasicFlow, we would be happy to include them in this section. Please open an issue for more arrangements or send a pull request with the benchmarks results.

benchmarks/rotatingDrum/readme.md

@@ -4,6 +4,12 @@
 
 This benchmark compares the performance of phasicFlow with a well-stablished commercial DEM software for simulating a rotating drum with varying particle counts (250k to 8M particles). The benchmark measures both computational efficiency and memory usage across different hardware configurations.
 
+**Summary of Results:**
+
+- phasicFlow achieves approximately 20% faster calculation than the commercial DEM software on the same hardware.
+- phasicFlow shows a 30% performance improvement when using the NVIDIA RTX A4000 compared to the RTX 4050Ti.
+- Memory usage is approximately 42% lower in phasicFlow compared to the commercial software, with phasicFlow using about 0.7 GB of memory per million particles, while the commercial software uses about 1.2 GB per million particles
+
 ## Simulation Setup
 
 <div align="center">

src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp

@@ -251,6 +251,8 @@ struct pwInteractionFunctor
 		realx3 xi = pos_[i];
 		
 		realx3x3 tri = triangles_(tj);
+		const realx3& normW = triangles_.normal(tj);
+
 		real ovrlp;
 		realx3 Nij, cp;
 
@@ -262,7 +264,7 @@ struct pwInteractionFunctor
 
 			int32 mInd	= wTriMotionIndex_[tj];
 			
-			auto Vw = motionModel_(mInd, cp);
+			auto Vw = motionModel_(mInd, cp, normW);
 			
 			//output<< "par-wall index "<< i<<" - "<< tj<<endl;
 			

src/Interaction/sphereInteraction/sphereInteraction/sphereInteractionKernels.hpp

@@ -238,7 +238,9 @@ struct pwInteractionFunctor
 		real Rj = 10000.0;
 		realx3 xi = pos_[i];
 		
-		realx3x3 tri = triangles_(tj); 
+		const realx3x3 tri = triangles_(tj);
+		const realx3& normW = triangles_.normal(tj); 
+
 		real ovrlp;
 		realx3 Nij, cp;
 
@@ -250,7 +252,7 @@ struct pwInteractionFunctor
 
 			int32 mInd	= wTriMotionIndex_[tj];
 			
-			auto Vw = motionModel_(mInd, cp);
+			auto Vw = motionModel_(mInd, cp, normW);
 			
 			//output<< "par-wall index "<< i<<" - "<< tj<<endl;
 			

src/MotionModel/MotionModel/MotionModel.hpp

@@ -85,15 +85,15 @@ public:
 		~ModelInterface()=default;
 
 		INLINE_FUNCTION_HD
-		realx3 pointVelocity(uint32 n, const realx3& p)const
+		realx3 pointVelocity(uint32 n, const realx3& p, const realx3& wallNormal)const
 		{
-			return components_[n].linVelocityPoint(p);
+			return components_[n].linVelocityPoint(p, wallNormal);
 		}
 
 		INLINE_FUNCTION_HD
-		realx3 operator()(uint32 n, const realx3& p)const
+		realx3 operator()(uint32 n, const realx3& p, const realx3& wallNormal)const
 		{
-			return pointVelocity(n,p);
+			return pointVelocity(n, p, wallNormal);
 		}
 
 		INLINE_FUNCTION_HD

src/MotionModel/entities/conveyorBelt/conveyorBelt.cpp

@@ -28,7 +28,7 @@ pFlow::conveyorBelt::conveyorBelt(const dictionary& dict)
   if(!read(dict))
 	{
 		fatalErrorInFunction<<
-		"  error in reading conveyorBelt from dictionary "<< dict.globalName()<<endl;
+		"  error in reading from dictionary "<< dict.globalName()<<endl;
 		fatalExit;
 	}
 }
@@ -37,7 +37,21 @@ FUNCTION_H
 bool pFlow::conveyorBelt::read(const dictionary& dict)
 {	
 
-	tangentVelocity_ 		= dict.getVal<realx3>("tangentVelocity");
+	linearVelocity_ 		= dict.getVal<real>("linearVelocity");
+    normal_             = dict.getVal<realx3>("normal");
+
+    if(normal_.length() > verySmallValue)
+    {
+        normal_.normalize();
+    }
+    else
+    {
+        fatalErrorInFunction<<
+            "  normal vector in "<< 
+            dict.globalName() << 
+            " cannot be zero vector "<<endl;
+        return false;
+    }
 	
 	return true;
 }
@@ -45,12 +59,19 @@ bool pFlow::conveyorBelt::read(const dictionary& dict)
 FUNCTION_H
 bool pFlow::conveyorBelt::write(dictionary& dict) const
 {
-	if( !dict.add("tangentVelocity", tangentVelocity_) )
+	if( !dict.add("linearVelocity", linearVelocity_) )
 	{
 		fatalErrorInFunction<<
 		"  error in writing tangentVelocity to dictionary "<< dict.globalName()<<endl;
 		return false;
 	}
+
+    if(!dict.add("normal", normal_))
+    {
+        fatalErrorInFunction<<
+        "  error in writing normal to dictionary "<< dict.globalName()<<endl;
+        return false;
+    }
 	return true;
 }
 
@@ -65,6 +86,7 @@ bool pFlow::conveyorBelt::read(iIstream& is)
 FUNCTION_H
 bool pFlow::conveyorBelt::write(iOstream& os)const
 {
-  os.writeWordEntry("tangentVelocity", tangentVelocity_);
-  return true;
+    os.writeWordEntry("linearVelocity", linearVelocity_);
+    os.writeWordEntry("normal", normal_);
+    return true;
 }

src/MotionModel/entities/conveyorBelt/conveyorBelt.hpp

@@ -41,65 +41,76 @@ class conveyorBelt
 {
 private:
 
-	realx3  tangentVelocity_{0, 0, 0};
+    /// @brief  linear velocity of the conveyor belt
+    real  linearVelocity_{0};
+
+    /// normal vector of the velocity plane.
+    /// The velocity vector is tangent to this plane (velocity plane). 
+    realx3 normal_{1,0,0};
 
 public:
 
-	TypeInfoNV("conveyorBelt");
+    TypeInfoNV("conveyorBelt");
 
-	FUNCTION_HD
-	conveyorBelt()=default;
+    FUNCTION_HD
+    conveyorBelt()=default;
 
-	FUNCTION_H
-	explicit conveyorBelt(const dictionary& dict);
+    FUNCTION_H
+    explicit conveyorBelt(const dictionary& dict);
     
 
-	FUNCTION_HD
-	conveyorBelt(const conveyorBelt&) = default;
+    FUNCTION_HD
+    conveyorBelt(const conveyorBelt&) = default;
 
-	conveyorBelt& operator=(const conveyorBelt&) = default;
+    conveyorBelt& operator=(const conveyorBelt&) = default;
 
-	INLINE_FUNCTION_HD
-	void setTime(real t)
-	{}
+    INLINE_FUNCTION_HD
+    void setTime(real t)
+    {}
 
-	INLINE_FUNCTION_HD
-	realx3 linVelocityPoint(const realx3 &)const
-	{
-		return tangentVelocity_;
-	}
+    /*INLINE_FUNCTION_HD
+    realx3 linVelocityPoint(const realx3 &)const
+    {
+        return tangentVelocity_;
+    }*/
 
-	INLINE_FUNCTION_HD
-	realx3 transferPoint(const realx3& p, real)const
-	{
-		return p;
-	}
+    INLINE_FUNCTION_HD
+    realx3 linVelocityPoint(const realx3 &, const realx3& wallNormal)const
+    {
+        return linearVelocity_ * cross(wallNormal, normal_);
+    }
 
-	// - IO operation
-	FUNCTION_H 
-	bool read(const dictionary& dict);
+    INLINE_FUNCTION_HD
+    realx3 transferPoint(const realx3& p, real)const
+    {
+        return p;
+    }
 
-	FUNCTION_H
-	bool write(dictionary& dict) const;
+    // - IO operation
+    FUNCTION_H 
+    bool read(const dictionary& dict);
 
-	FUNCTION_H
-	bool read(iIstream& is);
+    FUNCTION_H
+    bool write(dictionary& dict) const;
 
-	FUNCTION_H
-	bool write(iOstream& os)const;
+    FUNCTION_H
+    bool read(iIstream& is);
+
+    FUNCTION_H
+    bool write(iOstream& os)const;
 
 };
 
 inline iOstream& operator <<(iOstream& os, const conveyorBelt& obj)
 {
     
-	return os;
+    return os;
 }
 
 inline iIstream& operator >>(iIstream& is, conveyorBelt& obj)
 {
     
-	return is;
+    return is;
 }
 
 }

src/MotionModel/entities/multiRotatingAxis/multiRotatingAxis.hpp

@@ -120,7 +120,7 @@ public:
 
 		/// Tangential velocity at point p
 		INLINE_FUNCTION_HD
-		realx3 pointTangentialVel(const realx3& p)const
+		realx3 pointTangentialVel(const realx3& p, const realx3& wallNormal)const
 		{
 			realx3 parentVel(0);
 			auto parIndex = parentAxisIndex();
@@ -128,11 +128,11 @@ public:
 			while(parIndex != -1)
 			{
 				auto& ax = axisList_[parIndex];
-				parentVel += ax.linVelocityPoint(p);
+				parentVel += ax.linVelocityPoint(p, wallNormal);
 				parIndex = ax.parentAxisIndex();
 			}
 
-			return parentVel + rotatingAxis::linVelocityPoint(p);
+			return parentVel + rotatingAxis::linVelocityPoint(p, wallNormal);
 		}
 
 		/// Translate point p for dt seconds based on the axis information

src/MotionModel/entities/rotatingAxis/rotatingAxis.hpp

@@ -126,7 +126,7 @@ public:
 
 	/// Linear tangential velocity at point p
 	INLINE_FUNCTION_HD
-	realx3 linVelocityPoint(const realx3 &p)const;
+	realx3 linVelocityPoint(const realx3 &p, const realx3& wallNormal)const;
 
 	INLINE_FUNCTION_HD
 	realx3 transferPoint(const realx3 p, real dt)const;

src/MotionModel/entities/rotatingAxis/rotatingAxisI.hpp

@@ -20,7 +20,7 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 INLINE_FUNCTION_HD
-pFlow::realx3 pFlow::rotatingAxis::linVelocityPoint(const realx3 &p)const
+pFlow::realx3 pFlow::rotatingAxis::linVelocityPoint(const realx3 &p, const realx3&)const
 {
 	
 	if(!inTimeRange()) return {0,0,0};

src/MotionModel/entities/stationary/stationary.hpp

@@ -60,7 +60,7 @@ public:
 	{}
 
 	INLINE_FUNCTION_HD
-	realx3 linVelocityPoint(const realx3 &)const
+	realx3 linVelocityPoint(const realx3 &, const realx3&)const
 	{
 		return realx3(0);
 	}

src/MotionModel/entities/vibrating/vibrating.hpp

@@ -39,17 +39,17 @@ class dictionary;
  * Creates a sinoidal virating model for a wall. The viration is defined by
  * frequency, amplitude and phase angle. 
  * \f[ 
-	\vec{v}(t) = \vec{A} sin(\vec{\omega}(t-startTime)+\vec{\phi}) 
+    \vec{v}(t) = \vec{A} sin(\vec{\omega}(t-startTime)+\vec{\phi}) 
   \f]
  \verbatim
  // This creates sinoidal vibration on the wall in x-direction. The viration
  // starts at t = 0 s and ends at t = 10 s. 
  {
-  	angularFreq 	(10 0 0);
-  	amplitude 		( 1 0 0);
-  	phaseAngle		( 0 0 0);
-  	startTime 		0;
-  	endTime 		10;
+      angularFreq 	(10 0 0);
+      amplitude 		( 1 0 0);
+      phaseAngle		( 0 0 0);
+      startTime 		0;
+      endTime 		10;
  } \endverbatim
  * 
  * * 
@@ -64,102 +64,102 @@ class dictionary;
  */
 class vibrating
 :
-	public timeInterval
+    public timeInterval
 {
 
 private:
     
-	// rotation speed 
-	realx3	angularFreq_{0,0,0};
+    // rotation speed 
+    realx3	angularFreq_{0,0,0};
 
-	realx3 	phaseAngle_{0,0,0};
+    realx3 	phaseAngle_{0,0,0};
 
-	realx3  amplitude_{0,0,0};
+    realx3  amplitude_{0,0,0};
 
-	realx3  velocity_{0,0,0};
+    realx3  velocity_{0,0,0};
 
-	realx3  velocity0_{0,0,0};
+    realx3  velocity0_{0,0,0};
 
-	INLINE_FUNCTION_HD
-	void calculateVelocity()
-	{
-		if(inTimeRange())
-		{
-			velocity_ = amplitude_ * sin( angularFreq_*(time()-startTime() ) + phaseAngle_ );
-		}else
-		{
-			velocity_ = {0,0,0};
-		}
-	}
+    INLINE_FUNCTION_HD
+    void calculateVelocity()
+    {
+        if(inTimeRange())
+        {
+            velocity_ = amplitude_ * sin( angularFreq_*(time()-startTime() ) + phaseAngle_ );
+        }else
+        {
+            velocity_ = {0,0,0};
+        }
+    }
 
 public:
 
-	TypeInfoNV("vibrating");
+    TypeInfoNV("vibrating");
 
-	FUNCTION_HD
-	vibrating()=default;
+    FUNCTION_HD
+    vibrating()=default;
 
-	FUNCTION_H
-	explicit vibrating(const dictionary& dict);
+    FUNCTION_H
+    explicit vibrating(const dictionary& dict);
     
 
-	FUNCTION_HD
-	vibrating(const vibrating&) = default;
+    FUNCTION_HD
+    vibrating(const vibrating&) = default;
 
-	vibrating& operator=(const vibrating&) = default;
+    vibrating& operator=(const vibrating&) = default;
 
-	INLINE_FUNCTION_HD
-	void setTime(real t)
-	{
-		if( !equal(t, time()) ) velocity0_ = velocity_;
-		timeInterval::setTime(t);
-		calculateVelocity();
-	}
+    INLINE_FUNCTION_HD
+    void setTime(real t)
+    {
+        if( !equal(t, time()) ) velocity0_ = velocity_;
+        timeInterval::setTime(t);
+        calculateVelocity();
+    }
 
-	INLINE_FUNCTION_HD
-	realx3 linVelocityPoint(const realx3 &)const
-	{
-		return velocity_;
-	}
+    INLINE_FUNCTION_HD
+    realx3 linVelocityPoint(const realx3 &, const realx3&)const
+    {
+        return velocity_;
+    }
 
-	INLINE_FUNCTION_HD
-	realx3 transferPoint(const realx3& p, real dt)const
-	{
-		if(!inTimeRange()) return p;
-		return p + static_cast<real>(0.5*dt)*(velocity0_+velocity_);
-	}
+    INLINE_FUNCTION_HD
+    realx3 transferPoint(const realx3& p, real dt)const
+    {
+        if(!inTimeRange()) return p;
+        return p + static_cast<real>(0.5*dt)*(velocity0_+velocity_);
+    }
 
-	// - IO operation
-	FUNCTION_H 
-	bool read(const dictionary& dict);
+    // - IO operation
+    FUNCTION_H 
+    bool read(const dictionary& dict);
 
-	FUNCTION_H
-	bool write(dictionary& dict) const;
+    FUNCTION_H
+    bool write(dictionary& dict) const;
 
-	FUNCTION_H
-	bool read(iIstream& is);
+    FUNCTION_H
+    bool read(iIstream& is);
 
-	FUNCTION_H
-	bool write(iOstream& os)const;
+    FUNCTION_H
+    bool write(iOstream& os)const;
 
 };
 
 inline iOstream& operator <<(iOstream& os, const vibrating& obj)
 {
-	if(!obj.write(os))
-	{
-		fatalExit;
-	}
-	return os;
+    if(!obj.write(os))
+    {
+        fatalExit;
+    }
+    return os;
 }
 
 inline iIstream& operator >>(iIstream& is, vibrating& obj)
 {
-	if( !obj.read(is) )
-	{
-		fatalExit;
-	}
-	return is;
+    if( !obj.read(is) )
+    {
+        fatalExit;
+    }
+    return is;
 }
 
 }

src/PostprocessData/CMakeLists.txt

@@ -25,6 +25,8 @@ set(SourceFiles
     operation/PostprocessOperation/PostprocessOperationSum.cpp
     operation/PostprocessOperation/PostprocessOperationAverage.cpp
     operation/PostprocessOperation/PostprocessOperationAvMassVelocity.cpp
+    operation/PostprocessOperation/PostprocessOperationSolidVolFraction.cpp
+    operation/PostprocessOperation/PostprocessOperationBulkDensity.cpp 
        
     operation/includeMask/includeMask.cpp
     operation/includeMask/IncludeMasks.cpp

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAvMassVelocity.hpp

@@ -23,106 +23,10 @@ Licence:
 
 /*!
  * @class PostprocessOperationAvMassVelocity
- * @brief A class for averaging field values within specified regions during post-processing.
+ * @brief Calculates mass-weighted average velocity of particles in the regions
  * 
- * @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
+ * @see PostprocessOperationAverage
  */
 
 #include <variant>

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationAverage.cpp

@@ -162,6 +162,14 @@ bool PostprocessOperationAverage::write(const fileSystem &parDir) const
             processedFieldName()+"_prime2" + ".Start_" + ti.timeName());
         os2Ptr_ = makeUnique<oFstream>(path);
 
+        if(regPoints().scientific())
+        {
+            // set output format to scientific notation
+            os2Ptr_().stdStream()<<std::scientific;
+        }
+    
+        os2Ptr_().precision(regPoints().precision());
+
         regPoints().write(os2Ptr_());
     }
 

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationBulkDensity.cpp

@@ -0,0 +1,25 @@
+#include "PostprocessOperationBulkDensity.hpp"
+
+namespace pFlow::postprocessData
+{
+
+PostprocessOperationBulkDensity::PostprocessOperationBulkDensity
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    PostprocessOperationSum
+    (
+        opDict,
+        "mass",
+        "one",
+        "all",
+        regPoints,
+        fieldsDB
+    )
+{
+}
+
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationBulkDensity.hpp

@@ -0,0 +1,77 @@
+/*------------------------------- 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 __PostprocessOperationBulkDensity_hpp__
+#define __PostprocessOperationBulkDensity_hpp__
+
+/*!
+ * @class PostprocessOperationBulkDensity
+ * @brief Calculates bulk density in the regions
+ * 
+ *
+ * @see PostprocessOperationSum
+ */
+
+#include "PostprocessOperationSum.hpp"
+
+namespace pFlow::postprocessData
+{
+
+class PostprocessOperationBulkDensity
+:
+    public PostprocessOperationSum
+{
+
+public:
+
+    TypeInfo("PostprocessOperation<bulkDensity>");
+
+    /// @brief Constructs average operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
+    PostprocessOperationBulkDensity(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    /// destructor 
+    ~PostprocessOperationBulkDensity() override = default;
+
+    /// add this virtual constructor to the base class
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationBulkDensity,
+        dictionary
+    );
+
+    bool divideByVolume()const override
+    {
+        return true;
+    }
+
+};
+
+} 
+
+#endif //__PostprocessOperationSolidVolFraction_hpp__
+
+

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSolidVolFraction.cpp

@@ -0,0 +1,25 @@
+#include "PostprocessOperationSolidVolFraction.hpp"
+
+namespace pFlow::postprocessData
+{
+
+PostprocessOperationSolidVolFraction::PostprocessOperationSolidVolFraction
+(
+    const dictionary &opDict, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    PostprocessOperationSum
+    (
+        opDict,
+        "volume",
+        "one",
+        "all",
+        regPoints,
+        fieldsDB
+    )
+{
+}
+
+}

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSolidVolFraction.hpp

@@ -0,0 +1,77 @@
+/*------------------------------- 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 __PostprocessOperationSolidVolFraction_hpp__
+#define __PostprocessOperationSolidVolFraction_hpp__
+
+/*!
+ * @class PostprocessOperationSolidVolFraction
+ * @brief Calculates solid volume fraction in the regions
+ * 
+ *
+ * @see PostprocessOperationSum
+ */
+
+#include "PostprocessOperationSum.hpp"
+
+namespace pFlow::postprocessData
+{
+
+class PostprocessOperationSolidVolFraction
+:
+    public PostprocessOperationSum
+{
+
+public:
+
+    TypeInfo("PostprocessOperation<solidVolFraction>");
+
+    /// @brief Constructs average operation processor
+    /// @param opDict Operation parameters dictionary
+    /// @param regPoints Region points data
+    /// @param fieldsDB Fields database
+    PostprocessOperationSolidVolFraction(
+        const dictionary& opDict,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+    
+    /// destructor 
+    ~PostprocessOperationSolidVolFraction() override = default;
+
+    /// add this virtual constructor to the base class
+    add_vCtor
+    (
+        postprocessOperation,
+        PostprocessOperationSolidVolFraction,
+        dictionary
+    );
+
+    bool divideByVolume()const override
+    {
+        return true;
+    }
+
+};
+
+} 
+
+#endif //__PostprocessOperationSolidVolFraction_hpp__
+
+

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.cpp

@@ -41,6 +41,49 @@ PostprocessOperationSum::PostprocessOperationSum
     }
 }
 
+PostprocessOperationSum::PostprocessOperationSum
+(
+    const dictionary &opDict, 
+    const word &fieldName, 
+    const word &phiName, 
+    const word &includeName, 
+    const regionPoints &regPoints, 
+    fieldsDataBase &fieldsDB
+)
+:
+    postprocessOperation(
+        opDict, 
+        fieldName, 
+        phiName, 
+        includeName, 
+        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 PostprocessOperationSum::execute
 (

src/PostprocessData/operation/PostprocessOperation/PostprocessOperationSum.hpp

@@ -154,6 +154,14 @@ public:
         const regionPoints& regPoints,
         fieldsDataBase& fieldsDB);
     
+    PostprocessOperationSum(
+        const dictionary& opDict,
+        const word& fieldName,
+        const word& phiName,
+        const word& includeName,
+        const regionPoints& regPoints,
+        fieldsDataBase& fieldsDB);
+
     /// destructor 
     ~PostprocessOperationSum() override = default;
 

src/PostprocessData/operation/postprocessOperation/postprocessOperation.cpp

@@ -107,6 +107,14 @@ bool postprocessOperation::write(const fileSystem &parDir) const
             processedFieldName() + ".Start_" + ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
         
+        if(regPoints().scientific())
+        {   
+            // set output format to scientific notation 
+            osPtr_().stdStream()<<std::scientific;
+        }
+
+        osPtr_().precision(regPoints().precision());
+        
         regPoints().write(osPtr_());
     }
 

src/PostprocessData/operation/postprocessOperation/postprocessOperation.hpp

@@ -225,6 +225,7 @@ public:
     }
 
     /// whether the result is divided by volume of the region
+    virtual  
     bool divideByVolume()const
     {
         return divideByVolume_();

src/PostprocessData/postprocessComponent/PostprocessComponent/PostprocessComponent.cpp

@@ -157,6 +157,14 @@ bool pFlow::postprocessData::PostprocessComponent<RegionType, ProcessMethodType>
 
         auto osPtr = makeUnique<oFstream>(file);
 
+        // set output format to scientific notation
+        if(regPoints().scientific())
+        {
+            osPtr->stdStream() << std::scientific;
+        }
+    
+        osPtr().precision(regPoints().precision());
+
         regPoints().write(osPtr());
 
         for(auto& operation:operatios_)

src/PostprocessData/postprocessComponent/particleProbePostprocessComponent/particleProbePostprocessComponent.cpp

@@ -157,6 +157,12 @@ bool pFlow::postprocessData::particleProbePostprocessComponent::write(const file
         // file is not open yet
         fileSystem path = parDir + (name_+".Start_"+ti.timeName());
         osPtr_ = makeUnique<oFstream>(path);
+        
+        if(regionPointsPtr_().scientific())
+        {
+            osPtr_().stdStream() << std::scientific;
+        }
+        osPtr_().precision(regionPointsPtr_().precision());
         regionPointsPtr_().write(osPtr_());
     }
 

src/PostprocessData/readme.md

@@ -138,6 +138,41 @@ Regions define where in the domain the postprocessing operations are applied:
 | <td colspan="3">\* Particles selection is done when simulation reaches the time that is specified by `startTime` of the post-process component and this selection remains intact up to the end of simulation. This is very good for particle tracking purposes or when you want to analyze specific particles behavior over time.</td> |
 | <td colspan="3">\** This region creates a rectangular mesh and particles are located into cells according to their center points. When using `GaussianDistribution` as `processMethod`, a larger neighbor radius is considered for each cell and particles inside this neighbor radius are included in the calculations.</td> |
 
+### output format
+
+The output format of the postprocessing results can be controlled by the `precision` and `scientific` parameters:
+
+- `precision`: Number of decimal places for the output (defualt is 6).
+- `scientific`: Whether to use scientific notation for large numbers (options: `yes`, `no`, default is `yes`).
+
+for example, if you want to use 5 decimal places and no scientific notation, you can set:
+
+```C++
+on_single_sphere
+{
+    processMethod    arithmetic;
+
+    processRegion    sphere;
+    
+    sphereInfo
+    {
+        radius  0.01;
+        center  (-0.08 -0.08 0.015);
+    }
+    
+    timeControl default;
+    
+    precision    5; // default is 6
+
+    scientific   no; // default is yes
+
+    operations
+    (
+        // a list of operations should be defined here
+    );
+}
+```
+
 ## 6. Processing Operations for Bulk Properties
 
 Within each processing region of type `bulk`, you can define multiple operations to be performed:
@@ -185,6 +220,8 @@ In addition to the above basic functions, some derived functions are available f
 | Function | Property type | Description | Formula | Required Parameters |
 |----------|---------------|-------------|---------|---------------------|
 |`avMassVelocity` | bulk | Average velocity weighted by mass | $\frac{\sum_{i \in \text{region}} w_i \cdot m_i \cdot v_i}{\sum_{i \in \text{region}} w_i \cdot m_i}$ | - |
+|`solidVolFraction`| bulk| Volume fraction of solid| $\phi = \frac{\sum_{i \in \text{region}} w_i \cdot V_i}{V_{\text{region}}}$ | - |
+|`bulkDensity`| bulk| Bulk density of particles in the region | $\rho_{bulk} = \frac{\sum_{i \in \text{region}} w_i \cdot m_i}{V_{\text{region}}}$ | - |
 
 ### 6.4. Available Fields
 
@@ -566,6 +603,7 @@ components
         // are selected. Selection occurs at startTime: particles 
         // that are inside the box at t = startTime. 
         selector           box;
+
         boxInfo
         {
             min            (0 0 0);
@@ -575,6 +613,14 @@ components
         // center position of selected particles are processed 
         field              position;
 
+        // precision for output to file (optional: default is 6)
+        precision           8;
+
+        // if the output format of numbers in scientific format 
+        // is required, set scientific to yes, otherwise no
+        // (optional: default is yes)
+        scientific          no;
+        
         timeControl        simulationTime;
         // execution starts at 1.0 s
         startTime           1.0;
@@ -608,6 +654,14 @@ components
             cellExtension 3; 
         }
 
+        // precision for output to file (optional: default is 6)
+        precision   5;
+
+        // if the output format of numbers in scientific format
+        // is required, set scientific to yes, otherwise no
+        // (optional: default is yes)
+        scientific  no; 
+
         operations
         (
             avVelocity

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

@@ -12,7 +12,10 @@ regionPoints::regionPoints
 )
 :
     fieldsDataBase_(fieldsDataBase)
-{}
+{
+    precision_ = dict.getValOrSet<int>("precision", 6);
+    scientific_ = dict.getValOrSet<Logical>("scientific", Logical(true));
+}
 
 const Time& regionPoints::time() const
 {

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

@@ -54,6 +54,12 @@ class regionPoints
     /// Reference to the fields database containing simulation data
     fieldsDataBase& fieldsDataBase_;
 
+    /// default precision for output
+    int precision_ = 6;
+
+    /// if scientific notation is used for output
+    Logical scientific_;
+
 public:
 
     TypeInfo("regionPoints");
@@ -75,6 +81,15 @@ public:
     /// Returns non-const reference to the fields database
     fieldsDataBase& database();
 
+    int precision() const
+    {
+        return precision_;
+    }
+    
+    bool scientific()const
+    {
+        return scientific_(); 
+    }
     /// @brief size of elements 
     virtual 
     uint32 size()const = 0;

src/phasicFlow/triSurface/triSurface.hpp

@@ -46,6 +46,9 @@ private:
 	deviceViewType1D<realx3> dPoints_;
 	
 	deviceViewType1D<uint32x3>  dVectices_;
+
+	deviceViewType1D<realx3> dNormals_;
+
 public:
 
 	INLINE_FUNCTION_H
@@ -53,12 +56,14 @@ public:
 		uint32 	numPoints,
 		deviceViewType1D<realx3> points,
 		uint32 	numTrianlges,
-		deviceViewType1D<uint32x3> vertices )
+		deviceViewType1D<uint32x3> vertices,
+		deviceViewType1D<realx3> normals )
 	:
 		numPoints_(numPoints),
 		numTriangles_(numTrianlges),
 		dPoints_(points),
-		dVectices_(vertices)
+		dVectices_(vertices),
+		dNormals_(normals)
 	{}
 
 	INLINE_FUNCTION_HD
@@ -91,11 +96,17 @@ public:
 	INLINE_FUNCTION_HD
 	realx3x3 operator[](uint32 i)const { return triangle(i);	}
 
+	INLINE_FUNCTION_HD
+	const realx3& normal(uint32 i)const
+	{
+		return dNormals_[i];
+	}
+
 	INLINE_FUNCTION_HD
 	uint32 numPoints()const { return numPoints_; }
 
 	INLINE_FUNCTION_HD
-	uint32 numTrianlges()const { return numTriangles_;}
+	uint32 numTriangles()const { return numTriangles_;}
 };
 
 
@@ -224,7 +235,8 @@ public:
 				points_.size(),
 				points_.deviceViewAll(),
 				vertices_.size(),
-				vertices_.deviceViewAll());
+				vertices_.deviceViewAll(),
+				normals_.deviceViewAll() );
 		}
 		
 	//// - IO operations 

tutorials/sphereGranFlow/conveyorBelt/settings/geometryDict

@@ -7,14 +7,18 @@ objectType  dictionary;
 fileFormat  ASCII;
 /*---------------------------------------------------------------------------*/
 
-// motion model can be rotatingAxis or stationary or vibrating
 motionModel          conveyorBelt;               
 
 conveyorBeltInfo
 {
     conveyorBelt1
     {
-        tangentVelocity (0.5 0 0);
+        // linear velocity of belt
+        linearVelocity  0.5;
+
+        // normal vector of velocity plate
+        // The velocity vector is tangent to this plane  
+        normal          (0 -1 0);
     }
 }
 

tutorials/sphereGranFlow/conveyorBelt/settings/settingsDict

@@ -34,7 +34,7 @@ writeFormat                       ascii;          // data writting format (ascii
 
 timersReport                        Yes;          // report timers
 
-timersReportInterval               0.01;          // time interval for reporting timers
+timersReportInterval               0.1;          // time interval for reporting timers
 
 
 

tutorials/sphereGranFlow/screwConveyor/settings/geometryDict

@@ -38,7 +38,7 @@ surfaces
 		type 	 stlWall;  	// type of the wall
 		file 	 shell.stl;	// file name in stl folder		
 		material prop1;     // material name of this wall
-		motion 	 none;		// this surface is not moving ==> none
+		motion 	 none;		// this surface is not movng ==> none
 	}
 }
 

utilities/pFlowToVTK/pointFieldToVTK.cpp

@@ -225,8 +225,9 @@ bool pFlow::PFtoVTK::addUndstrcuturedGridField(
 
     os << "DATASET UNSTRUCTURED_GRID\n";
     
-    if (numPoints == 0)
-        return true;
+    // this is commented to resolve the problem of no particle in paraview
+    //if (numPoints == 0)
+    //    return true;
 
     os << "POINTS " << numPoints << " float"<<'\n';
     if(bindaryOutput__)
@@ -349,8 +350,8 @@ bool pFlow::PFtoVTK::addRealPointField(
     const real *field,
     uint32 numData)
 {
-    if (numData == 0)
-        return true;
+    //if (numData == 0)
+    //    return true;
 
     os << "FIELD FieldData 1\n"
        << fieldName << " 1 " << numData << " float\n";
@@ -380,8 +381,8 @@ bool pFlow::PFtoVTK::addRealx3PointField(
     const realx3 *field,
     uint32 numData)
 {
-    if (numData == 0)
-        return true;
+    //if (numData == 0)
+    //    return true;
 
     os << "FIELD FieldData 1\n"
        << fieldName << " 3 " << numData << " float\n";

utilities/pFlowToVTK/pointFieldToVTK.hpp

@@ -131,8 +131,8 @@ namespace pFlow::PFtoVTK
 		uint32 numData)
 	{
 
-		if (numData == 0)
-			return true;
+		//if (numData == 0)
+		//	return true;
 
 		if(std::is_same_v<IntType, int> || std::is_same_v<IntType, const int> )
 		{

utilities/particlesPhasicFlow/positionFile/positionFile.cpp

@@ -49,7 +49,14 @@ bool pFlow::positionFile::positionPointsFile()
 		is >> tok;
 		if(tok.good()&& tok.isNumber()&& !is.eof())
 		{
-			tempPoint.x() = tok.realToken();
+			if(tok.isReal())
+			{
+				tempPoint.x() = tok.realToken();
+			}
+			else
+			{
+				tempPoint.x() = static_cast<real>(tok.int64Token());
+			}
 		}
 		else
 		{
@@ -71,7 +78,14 @@ bool pFlow::positionFile::positionPointsFile()
 		is >> tok;
 		if(tok.good()&& tok.isNumber()&& !is.eof())
 		{
-			tempPoint.y() = tok.realToken();
+			if(tok.isReal())
+			{
+				tempPoint.y() = tok.realToken();
+			}
+			else
+			{
+				tempPoint.y() = static_cast<real>(tok.int64Token());
+			}
 		}
 		else
 		{
@@ -93,7 +107,14 @@ bool pFlow::positionFile::positionPointsFile()
 		is >> tok;
 		if(tok.good()&& tok.isNumber()&& !is.eof())
 		{
-			tempPoint.z() = tok.realToken();
+			if(tok.isReal())
+			{
+				tempPoint.z() = tok.realToken();
+			}
+			else
+			{
+				tempPoint.z() = static_cast<real>(tok.int64Token());
+			}
 		}
 		else
 		{