Merge branch 'main' into documentation

CMakeLists.txt

@@ -77,7 +77,7 @@ add_subdirectory(solvers)
 add_subdirectory(utilities)
 
 add_subdirectory(DEMSystems)
-#add_subdirectory(test_newFeatures)
+#add_subdirectory(testIO)
 
 
 install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H"

doc/header.html

@@ -4,7 +4,7 @@
 <meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
 <meta http-equiv="X-UA-Compatible" content="IE=9"/>
 <meta name="generator" content="Doxygen $doxygenversion"/>
-<meta name="description" content="PhasicFlow is an open-source parallel DEM package for simulating granular flow developed in C++ and can be exectued on both GPU (like cuda) and CPU.">
+<meta name="description" content="PhasicFlow is an open-source parallel DEM (discrete element method) package for simulating granular flow. It is developed in C++ and can be exectued on both GPU (like CUDA) and CPU.">
 <!--BEGIN PROJECT_NAME-->
 <title>$projectname: $title</title><!--END PROJECT_NAME-->
 <!--BEGIN !PROJECT_NAME-->

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -28,45 +28,68 @@ Licence:
 namespace pFlow
 {
 
-
+/**
+ * Second order Adams-Bashforth integration method for solving ODE
+ * 
+ * This is a one-step integration method and does not have prediction step.
+ * 
+ */
 class AdamsBashforth2
 :
 	public integration
 {
 protected:
 
+	/// dy at t-dt
 	realx3PointField_D& dy1_;
 
+	/// Range policy for integration kernel (alias)
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
 			Kokkos::IndexType<int32>
 			>;
+
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsBashforth2");
 
-	//// - Constructors
+	// - Constructors
+
+		/// Construct from components
 		AdamsBashforth2(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth2>(*this);
+		}
+
+		/// Destructor 
 		virtual ~AdamsBashforth2()=default;
 
-		// - add a virtual constructor 
+		/// Add this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth2,
 			word);
 
 
-	//// - Methods
+	// - Methods
-		bool predict(real UNUSED(dt), realx3Vector_D& UNUSED(y), realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+		bool predict(
+			real UNUSED(dt), 
+			realx3Vector_D& UNUSED(y), 
+			realx3Vector_D& UNUSED(dy)) override;
+
+		bool correct(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -77,15 +100,20 @@ public:
 			return false;
 		}
 		
-		uniquePtr<integration> clone()const override
+		/// Integrate on all points in the active range 
-		{
+		bool intAll(
-			return makeUnique<AdamsBashforth2>(*this);
+			real dt, 
-		}
+			realx3Vector_D& y, 
-
+			realx3Vector_D& dy, 
-		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+			range activeRng);
 
+		/// Integrate on active points in the active range 
 		template<typename activeFunctor>
-		bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP );
+		bool intRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP );
 
 };
 

src/Integration/AdamsBashforth3/AdamsBashforth3.cpp

@@ -31,26 +31,8 @@ pFlow::AdamsBashforth3::AdamsBashforth3
 )
 :
 	integration(baseName, owner, pStruct, method),
-	/*dy1_(
-		owner.emplaceObject<realx3PointField_D>(
-			objectFile(
-				groupNames(baseName,"dy1"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			zero3)),
-	dy2_(
-		owner.emplaceObject<realx3PointField_D>(
-			objectFile(
-				groupNames(baseName,"dy2"),
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS),
-			pStruct,
-			zero3))*/
 	history_(
-		owner.emplaceObject<HistoryFieldType>(
+		owner.emplaceObject<pointField<VectorSingle,AB3History>>(
 			objectFile(
 				groupNames(baseName,"AB3History"),
 				"",

src/Integration/AdamsBashforth3/AdamsBashforth3.hpp

@@ -30,10 +30,10 @@ namespace pFlow
 
 struct AB3History
 {
+	TypeInfoNV("AB3History");
+
 	realx3 dy1_={0,0,0};
 	realx3 dy2_={0,0,0};
-
-	TypeInfoNV("AB3History");
 };
 
 
@@ -65,21 +65,21 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
 	return str;
 }
 
-
+/**
+ * Third order Adams-Bashforth integration method for solving ODE
+ * 
+ * This is a one-step integration method and does not have prediction step.
+ */
 class AdamsBashforth3
 :
 	public integration
 {
 protected:
 	
-	using HistoryFieldType = pointField<VectorSingle,AB3History>;
+	/// Integration history
-	//realx3PointField_D& dy1_;
+	pointField<VectorSingle,AB3History>& history_;
-
-	//realx3PointField_D& dy2_;
-
-	//  this is a device 
-	HistoryFieldType& history_;
 
+	/// Range policy for integration kernel
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -91,23 +91,32 @@ public:
 	// type info
 	TypeInfo("AdamsBashforth3");
 
-	//// - Constructors
+	// - Constructors
+		
+		/// Construct from components
 		AdamsBashforth3(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 		
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth3>(*this);
+		}
+
+		/// Destructor 
 		virtual ~AdamsBashforth3()=default;
 
-		// - add a virtual constructor 
+		/// Add this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth3,
 			word);
 
 
-	//// - Methods
+	// - Methods
+
 		bool predict(
 			real UNUSED(dt),
 			realx3Vector_D & UNUSED(y),
@@ -126,18 +135,17 @@ public:
 			return false;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Integrate on all points in the active range 
-		{
+		bool intAll(
-			return makeUnique<AdamsBashforth3>(*this);
+			real dt, 
-		}
-
-		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
 			range activeRng);
 
+		/// Integrate on active points in the active range 
 		template<typename activeFunctor>
-		bool intRange(real dt,
+		bool intRange(
+			real dt,
 			realx3Vector_D& y,
 			realx3Vector_D& dy,
 			activeFunctor activeP );

src/Integration/AdamsBashforth4/AdamsBashforth4.cpp

@@ -32,7 +32,7 @@ pFlow::AdamsBashforth4::AdamsBashforth4
 :
 	integration(baseName, owner, pStruct, method),
 	history_(
-		owner.emplaceObject<HistoryFieldType>(
+		owner.emplaceObject<pointField<VectorSingle,AB4History>>(
 			objectFile(
 				groupNames(baseName,"AB4History"),
 				"",

src/Integration/AdamsBashforth4/AdamsBashforth4.hpp

@@ -30,11 +30,12 @@ namespace pFlow
 
 struct AB4History
 {
+	TypeInfoNV("AB4History");
+
 	realx3 dy1_={0,0,0};
 	realx3 dy2_={0,0,0};
 	realx3 dy3_={0,0,0};
 	
-	TypeInfoNV("AB4History");
 };
 
 
@@ -68,23 +69,21 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
 	return str;
 }
 
-
+/**
+ * Fourth order Adams-Bashforth integration method for solving ODE
+ * 
+ * This is a one-step integration method and does not have prediction step.
+ */
 class AdamsBashforth4
 :
 	public integration
 {
 protected:
 
-	using HistoryFieldType = pointField<VectorSingle,AB4History>;
+	/// Integration history
-	//realx3PointField_D& dy1_;
+	pointField<VectorSingle,AB4History>& history_;
-
-	//realx3PointField_D& dy2_;
-
-
-
-	//  this is a device 
-	HistoryFieldType& history_;
 
+	/// Range policy for integration kernel
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -93,32 +92,42 @@ protected:
 
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsBashforth4");
 
-	//// - Constructors
+	// - Constructors
+		
+		/// Construct from components
 		AdamsBashforth4(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 		
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth4>(*this);
+		}
+
+		/// Destructor
 		virtual ~AdamsBashforth4()=default;
 
-		// - add a virtual constructor 
+		/// Add a this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth4,
 			word);
 
 
-	//// - Methods
+	// - Methods
+		
 		bool predict(
 			real UNUSED(dt),
 			realx3Vector_D & UNUSED(y),
 			realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, 
+		bool correct(
+			real dt, 
 			realx3Vector_D & y,
 			realx3Vector_D& dy) override;
 
@@ -131,18 +140,17 @@ public:
 			return false;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Integrate on all points in the active range 
-		{
+		bool intAll(
-			return makeUnique<AdamsBashforth4>(*this);
+			real dt, 
-		}
-
-		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
 			range activeRng);
 
+		/// Integrate on active points in the active range 
 		template<typename activeFunctor>
-		bool intRange(real dt,
+		bool intRange(
+			real dt,
 			realx3Vector_D& y,
 			realx3Vector_D& dy,
 			activeFunctor activeP );

src/Integration/AdamsBashforth5/AdamsBashforth5.cpp

@@ -32,7 +32,7 @@ pFlow::AdamsBashforth5::AdamsBashforth5
 :
 	integration(baseName, owner, pStruct, method),
 	history_(
-		owner.emplaceObject<HistoryFieldType>(
+		owner.emplaceObject<pointField<VectorSingle,AB5History>>(
 			objectFile(
 				groupNames(baseName,"AB5History"),
 				"",

src/Integration/AdamsBashforth5/AdamsBashforth5.hpp

@@ -30,12 +30,12 @@ namespace pFlow
 
 struct AB5History
 {
+	TypeInfoNV("AB5History");
+
 	realx3 dy1_={0,0,0};
 	realx3 dy2_={0,0,0};
 	realx3 dy3_={0,0,0};
 	realx3 dy4_={0,0,0};
-
-	TypeInfoNV("AB5History");
 };
 
 
@@ -71,18 +71,21 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
 	return str;
 }
 
-
+/**
+ * Fifth order Adams-Bashforth integration method for solving ODE
+ * 
+ * This is a one-step integration method and does not have prediction step.
+ */
 class AdamsBashforth5
 :
 	public integration
 {
 protected:
 
-	using HistoryFieldType = pointField<VectorSingle,AB5History>;
+	/// Integration history
-	
+	pointField<VectorSingle,AB5History>& history_;
-	//  this is a device 
-	HistoryFieldType& history_;
 
+	/// Range policy for integration kernel
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -91,32 +94,40 @@ protected:
 
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsBashforth5");
 
-	//// - Constructors
+	// - Constructors
+		
 		AdamsBashforth5(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth5>(*this);
+		}
+
 		virtual ~AdamsBashforth5()=default;
 
-		// - add a virtual constructor 
+		/// Add this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsBashforth5,
 			word);
 
 
-	//// - Methods
+	// - Methods
+		
 		bool predict(
 			real UNUSED(dt),
 			realx3Vector_D & UNUSED(y),
 			realx3Vector_D& UNUSED(dy)) override;
 
-		bool correct(real dt, 
+		bool correct(
+			real dt, 
 			realx3Vector_D & y,
 			realx3Vector_D& dy) override;
 
@@ -129,18 +140,17 @@ public:
 			return false;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Integrate on all points in the active range 
-		{
+		bool intAll(
-			return makeUnique<AdamsBashforth5>(*this);
+			real dt, 
-		}
-
-		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
 			range activeRng);
 
+		/// Integrate on active points in the active range 
 		template<typename activeFunctor>
-		bool intRange(real dt,
+		bool intRange(
+			real dt,
 			realx3Vector_D& y,
 			realx3Vector_D& dy,
 			activeFunctor activeP );
@@ -184,4 +194,4 @@ bool pFlow::AdamsBashforth5::intRange(
 
 } // pFlow
 
-#endif //__integration_hpp__
+#endif //

src/Integration/AdamsMoulton3/AdamsMoulton3.hpp

@@ -28,19 +28,27 @@ Licence:
 namespace pFlow
 {
 
-
+/**
+ * Third order Adams-Moulton integration method for solving ODE
+ * 
+ * This is a predictor-corrector integration method.
+ */
 class AdamsMoulton3
 :
 	public integration
 {
 protected:
 
+	/// y at time t
 	realx3PointField_D& y0_;
 
+	/// dy at time t
 	realx3PointField_D& dy0_;
 
+	/// dy at time t-dt
 	realx3PointField_D& dy1_;
 
+	/// Range policy for integration kernel
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -48,29 +56,44 @@ protected:
 			>;
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsMoulton3");
 
-	//// - Constructors
+	// - Constructors
+
+		/// Construct from components
 		AdamsMoulton3(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsMoulton3>(*this);
+		}
+
+		/// Destructor 
 		virtual ~AdamsMoulton3()=default;
 
-		// - add a virtual constructor 
+		/// Add this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsMoulton3,
 			word);
 
 
-	//// - Methods
+	// - Methods
-		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
 
-		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+		bool predict(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
+
+		bool correct(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -81,20 +104,35 @@ public:
 			return true;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Prediction step on all points in the active range
-		{
+		bool predictAll(
-			return makeUnique<AdamsMoulton3>(*this);
+			real dt, 
-		}
+			realx3Vector_D& y, 
-
+			realx3Vector_D& dy, 
-		bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+			range activeRng);
 
+		/// Prediction step on active points in the active range
 		template<typename activeFunctor>
-		bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP);
+		bool predictRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP);
 
-		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+		/// Integrate on all points in the active range
+		bool intAll(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			range activeRng);
 
+		/// Integrate on active points in the active range
 		template<typename activeFunctor>
-		bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP );
+		bool intRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP);
 
 };
 
@@ -176,4 +214,4 @@ bool pFlow::AdamsMoulton3::intRange(
 
 } // pFlow
 
-#endif //__integration_hpp__
+#endif //

src/Integration/AdamsMoulton4/AdamsMoulton4.hpp

@@ -28,21 +28,30 @@ Licence:
 namespace pFlow
 {
 
-
+/**
+ * Fourth order Adams-Moulton integration method for solving ODE
+ * 
+ * This is a predictor-corrector integration method.
+ */
 class AdamsMoulton4
 :
 	public integration
 {
 protected:
 
+	/// y at time t
 	realx3PointField_D& y0_;
 
+	/// dy at time t
 	realx3PointField_D& dy0_;
 
+	/// dy at time t-dt
 	realx3PointField_D& dy1_;
 
+	/// dy at time t-2*dt
 	realx3PointField_D& dy2_;
 
+	/// Range policy for integration kernel 
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -50,29 +59,44 @@ protected:
 			>;
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsMoulton4");
 
-	//// - Constructors
+	// - Constructors
+		
+		/// Construct from components
 		AdamsMoulton4(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsMoulton4>(*this);
+		}
+
+		/// Destructor 
 		virtual ~AdamsMoulton4()=default;
 
-		// - add a virtual constructor 
+		// Add this to the virtual constructor table 
 		add_vCtor(
 			integration,
 			AdamsMoulton4,
 			word);
 
 
-	//// - Methods
+	// - Methods
-		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
 
-		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+		bool predict(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
+
+		bool correct(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -83,20 +107,35 @@ public:
 			return true;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Prediction step on all points in the active range
-		{
+		bool predictAll(
-			return makeUnique<AdamsMoulton4>(*this);
+			real dt, 
-		}
+			realx3Vector_D& y, 
-
+			realx3Vector_D& dy, 
-		bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+			range activeRng);
 
+		/// Prediction step on active points in the active range
 		template<typename activeFunctor>
-		bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP);
+		bool predictRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP);
 
-		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+		/// Integrate on all points in the active range
+		bool intAll(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			range activeRng);
 
+		/// Integrate on active points in the active range
 		template<typename activeFunctor>
-		bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP );
+		bool intRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP );
 
 };
 
@@ -182,4 +221,4 @@ bool pFlow::AdamsMoulton4::intRange(
 
 } // pFlow
 
-#endif //__integration_hpp__
+#endif //

src/Integration/AdamsMoulton5/AdamsMoulton5.hpp

@@ -28,23 +28,33 @@ Licence:
 namespace pFlow
 {
 
-
+/**
+ * Fifth order Adams-Moulton integration method for solving ODE
+ * 
+ * This is a predictor-corrector integration method.
+ */
 class AdamsMoulton5
 :
 	public integration
 {
 protected:
 
+	/// y at time t
 	realx3PointField_D& y0_;
 
+	/// dy at time t
 	realx3PointField_D& dy0_;
 
+	/// dy at time t-dt
 	realx3PointField_D& dy1_;
 
+	/// dy at time t-2*dt
 	realx3PointField_D& dy2_;
 
+	/// dy at time t-3*dt
 	realx3PointField_D& dy3_;
 
+	/// Range policy for integration kernel 
 	using rpIntegration = Kokkos::RangePolicy<
 			DefaultExecutionSpace,
 			Kokkos::Schedule<Kokkos::Static>,
@@ -52,16 +62,24 @@ protected:
 			>;
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("AdamsMoulton5");
 
-	//// - Constructors
+	// - Constructors
+		
+		/// Construct from components
 		AdamsMoulton5(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 		
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsMoulton5>(*this);
+		}
+
+		/// Destructor 
 		virtual ~AdamsMoulton5()=default;
 
 		// - add a virtual constructor 
@@ -71,10 +89,17 @@ public:
 			word);
 
 
-	//// - Methods
+	// - Methods
-		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
 		
-		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+		bool predict(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
+
+		bool correct(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy) override;
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,
@@ -85,20 +110,35 @@ public:
 			return true;
 		}
 
-		uniquePtr<integration> clone()const override
+		/// Prediction step on all points in the active range
-		{
+		bool predictAll(
-			return makeUnique<AdamsMoulton5>(*this);
+			real dt, 
-		}
+			realx3Vector_D& y, 
-
+			realx3Vector_D& dy, 
-		bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+			range activeRng);
 
+		/// Prediction step on active points in the active range
 		template<typename activeFunctor>
-		bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP);
+		bool predictRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP);
 
-		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+		/// Integrate on all points in the active range
+		bool intAll(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			range activeRng);
 
+		/// Integrate on active points in the active range
 		template<typename activeFunctor>
-		bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP );
+		bool intRange(
+			real dt, 
+			realx3Vector_D& y, 
+			realx3Vector_D& dy, 
+			activeFunctor activeP );
 
 };
 

src/Integration/integration/integration.hpp

@@ -31,32 +31,71 @@ Licence:
 namespace pFlow
 {
 
-
+/**
+ * Base class for integrating the first order ODE (IVP) 
+ * 
+ * The ODE should be in the following form:
+ *\f[ 
+	 \frac{dy}{dt} = f(y,t)
+  \f]
+ * for example the equation of motion is in the following form:
+ *\f[ 
+	m\frac{d\vec{v}}{dt} = m\vec{g} + \sum \vec{f_c}(\vec{v},t)
+  \f]
+ *
+ * The integration method can be either one-step or predictor-corrector type.
+ * 
+ */
 class integration
 {
 protected:
 
+	// - Protected data members 
+
+		/// The owner repository that all fields are storred in
 		repository& owner_;
 
+		/// The base name for integration 
 		const word baseName_;
 
+		/// A reference to pointStructure
 		const pointStructure& pStruct_;
 
 public:
 
-	// type info
+	/// Type info
 	TypeInfo("integration");
 
-	//// - Constructors
+
+	// - Constructors
+
+		/// Construct from components
 		integration(
 			const word& baseName,
 			repository& owner,
 			const pointStructure& pStruct,
 			const word& method);
 
+		/// Copy constructor 
+		integration(const integration&) = default;
+
+		/// Move constructor
+		integration(integration&&) = default;
+
+		/// Copy assignment 
+		integration& operator = (const integration&) = default;
+
+		/// Move assignment
+		integration& operator = (integration&&) = default;
+
+		/// Polymorphic copy/cloning
+		virtual 
+		uniquePtr<integration> clone()const=0;
+
+		/// Destructor 
 		virtual ~integration()=default;
 
-		// - add a virtual constructor 
+		/// Add a virtual constructor 
 		create_vCtor(
 			integration,
 			word,
@@ -67,35 +106,49 @@ public:
 			(baseName, owner, pStruct, method) );
 
 
-	//// - Methods
+	// - Methods
 
+		/// Const ref to pointStructure
+		inline
 		const auto& pStruct()const
 		{
 			return pStruct_;
 		}
 
-		virtual bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
+		/// Base name
-
+		inline
-		virtual bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
-
-		virtual bool setInitialVals(
-				const int32IndexContainer& newIndices,
-				const realx3Vector& y) = 0;
-
-		virtual bool needSetInitialVals()const = 0;
-
-		virtual uniquePtr<integration> clone()const=0;
-
 		const word& baseName()const
 		{
 			return baseName_;
 		}
 
+		/// Ref to the owner repository
+		inline
 		repository& owner()
 		{
 			return owner_;
 		}
 
+		/// Prediction step in integration
+		virtual 
+		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
+
+		/// Correction/main integration step
+		virtual 
+		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
+
+		/// Set the initial values for new indices 
+		virtual 
+		bool setInitialVals(
+				const int32IndexContainer& newIndices,
+				const realx3Vector& y) = 0;
+
+		/// Check if the method requires any set initial vals
+		virtual 
+		bool needSetInitialVals()const = 0;
+
+
+	/// Create the polymorphic object based on inputs
 	static
 		uniquePtr<integration> create(
 			const word& baseName,
@@ -103,7 +156,7 @@ public:
 			const pointStructure& pStruct,
 			const word& method);
 
-};
+}; // integration
 
 } // pFlow
 

src/Integration/integration/integrations.hpp

@@ -1,28 +0,0 @@
-/*------------------------------- 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 __integrations_hpp__
-#define __integrations_hpp__
-
-#include "integration.hpp"
-#include "AdamsBashforth2.hpp"
-#include "AdamsBashforth3.hpp"
-
-#endif

src/Particles/dynamicPointStructure/dynamicPointStructure.hpp

@@ -24,7 +24,7 @@ Licence:
 
 #include "Time.hpp"
 #include "pointFields.hpp"
-#include "integrations.hpp"
+#include "integration.hpp"
 #include "uniquePtr.hpp"
 
 namespace pFlow
@@ -57,17 +57,6 @@ public:
 
 	dynamicPointStructure(const dynamicPointStructure& ps) = default;
 
-	/*dynamicPointStructure(const dynamicPointStructure& ps):
-		pointStructure(ps),
-		time_(ps.time_),
-		integrationMethod_(ps.integrationMethod_),
-		velocity_(ps.velocity_),
-		integrationPos_(ps.integrationPos_->clone()),
-		integrationVel_(ps.integrationVel_->clone())
-	{
-
-	}*/
-
 		
 		// - no move construct 
 	dynamicPointStructure(dynamicPointStructure&&) = delete;

src/demComponent/demComponent.hpp

@@ -28,22 +28,37 @@ Licence:
 namespace pFlow
 {
 
+/**
+ * A base class for every main component of DEM system.
+ * 
+ * This class provides abstraction at a very high level for any solver/utility
+ * that forces the derived component (i.e. particles, geometry, and etc) to 
+ * advance over time when iterate is called in time a loop.
+ * 
+ */
 class demComponent
 {
 protected:
 
+	// - Protected data members 
+
+		/// Name of the DEM component 
 		word 	componentName_;
 
-
+		/// Reference to systemControl 
 		systemControl& control_;
 
-
+		/// All timers (if any) of this component
 		Timers 	timers_;
 
 public:
 
+	/// Type info
 	TypeInfo("demComponent");
 
+	// - Constructors
+	
+		/// construct from components
 		demComponent(const word& name, systemControl& control)
 		:
 			componentName_(name),
@@ -51,47 +66,94 @@ public:
 			timers_(name, &control.timers())
 		{}
 
+		/// No copy constructor
+		demComponent(const demComponent&) = delete;
+
+		/// No move constructor 
+		demComponent(demComponent&&) = delete;
+
+		/// No copy assignment
+		demComponent& operator = (const demComponent&) = delete;
+
+		/// No move assignment 
+		demComponent& operator =(demComponent&&) = delete;
+
+		/// destructor 
 		virtual ~demComponent() = default;
 
 	
+	// - Member functions 
+
+		/// Const ref to systemControl
+		inline
 		const auto& control()const
 		{
 			return control_;
 		}
 
+		/// Ref to systemControl
+		inline
 		auto& control()
 		{
 			return control_;
 		}
 
+		/// Time step of integration
 		inline
 		real dt()const
 		{
 			return control_.time().dt();
 		}
 
+		/// Current simulation time 
 		inline
 		real currentTime()const
 		{
 			return control_.time().currentTime();
 		}
 
-	auto& timers(){
+		/// Const ref to timers
+		inline
+		const auto& timers() const
+		{
 			return timers_;
 		}
 
-	const auto& timers() const{
+		/// Ref to timers
+		inline
+		auto& timers()
+		{
 			return timers_;
 		}
 
+		/// This is called before the start of time loop
+		virtual
+		bool beforeTimeLoop()
+		{
+			notImplementedFunction;
+			return false;
+		}
 
-	virtual bool beforeIteration() = 0;
+		/// This is called in time loop, before iterate
+		virtual 
+		bool beforeIteration() = 0;
 
+		/// This is called in time loop. Perform the main calculations 
+		/// when the component should evolve along time.
+		virtual 
+		bool iterate() = 0;
 
-	virtual bool iterate() = 0;
+		/// This is called in time loop, after iterate.
+		virtual 
+		bool afterIteration() = 0;
 
-
+		/// This is called after the time loop
-	virtual bool afterIteration() = 0;
+		virtual
+		bool afterTimeLoop()
+		{
+			notImplementedFunction;
+			return false;
+		}
 
 };
 

tutorials/sphereGranFlow/RotatingDrumWithBaffles/ReadMe.md

@@ -1,5 +1,5 @@
 # Problem Definition 
-The problem is to simulate a rotating drum with the diameter **0.24 m**, the length **0.1 m** and **6** Baffles, rotating at **15 rpm**. This drum is filled with **20000** Particles.The timestep for integration is **0.00001 s**. There are 2 types of Particles in this drum each are beining inserted during simulation to fill the drum.
+The problem is to simulate a rotating drum with the diameter **0.24 m**, the length **0.1 m** and **6** Baffles, rotating at **15 rpm**. This drum is filled with **20000** Particles.The timestep for integration is **0.00001 s**. There are 2 types of Particles in this drum each are being inserted during simulation to fill the drum.
 * **12500** Particles with **4 mm** diameter, at the rate of 12500 particles/s for 1 sec. 
 * **7500** Particles with **5mm** diameter, at the rate of 7500 particles/s for 1 sec.   
 

tutorials/sphereGranFlow/toteblender/ReadMe.md

@@ -1,6 +1,6 @@
 # Problem Definition
-The problem is to simulate a double pedestal tote blender with the diameter **0.03 m** and **0.1 m** respectively, the length **0.3 m**, rotating at **28 rpm**. This blender is filled with **20000** Particles. The timestep for integration is **0.00001 s**. There is one type of Particle in this blender that are being inserted during simulation to fill the blender.
+The problem is to simulate a double pedestal tote blender (mixer) with the diameter **0.03 m** and **0.1 m** respectively, the length **0.3 m**, rotating at **28 rpm**. This blender is filled with **24000** particles. The timestep for integration is **0.00001 s**. There is one type of particle in this blender. Particles are positioned before start of simulation to fill the blender.
-* **20000** particles with **4 mm** diameter, at the rate of 20000 particles/s for 1 sec. َAfter settling particles, this blender starts to rotate at t=**1s**. 
+* **24000** particles with **5 mm** diameter are positioned, in order, and let to be settled under gravity. After settling particles, this blender starts to rotate at t=**1s**. 
 
 <html>
 <body>
@@ -8,238 +8,22 @@ The problem is to simulate a double pedestal tote blender with the diameter **0.
 	a view of the tote-blender while rotating 
 </div></b>
 <div align="center">
-<img src="sample sample sample sample", width=700px>
+<img src="https://github.com/PhasicFlow/phasicFlow/blob/media/media/Tote-blender.gif", width=700px>
 </div>
+<div align="center"><i>
+	particles are colored according to their velocity
+</div></i>
 </body>
 </html>
 
 # Setting up the Case
-As it has been explained in the previous cases, the simulation case setup is based on text-based scripts. Here, the simulation case setup are sotred in two folders: `caseSetup`, `setting`. (see the above folders). Unlike the previous cases, this case does not have the `stl` file. and the geometry is described in the `geometryDict` file.
+As it has been explained in the previous cases, the simulation case setup is based on text-based scripts. Here, the simulation case setup files are stored into two folders: `caseSetup`, `setting` (see the above folders). Unlike the previous cases, this case does not have the `stl` file and the surfaces are defined based on the built-in utilities in phasicFlow. See next the section for more information on how we can setup the geometry and its rotation. 
 
-## Defining particles 
+## Geometry 
-Then in the `caseSetup/sphereShape` the diameter and the material name of the particles are defined.
-```C++
-// names of shapes 
-names 		(sphere1); 	
-// diameter of shapes (m)
-diameters 	(0.004);
-// material names for shapes 	
-materials	(prop1);	
-```
-## Particle Insertion
-In this case we have a region for ordering particles. These particles are placed in this blender. For example the script for the inserted particles is shown below.
 
-<div align="center"> 
+### Defining rotation axis
-in <b>caseSetup/particleInsertion</b> file
+In file `settings/geometryDict` the information of rotating axis and speed of rotation are defined. The rotation of this blender starts at time=**0.5 s** and ends at time=**9.5 s**.
-</div>
 
-```C++
-// positions particles 
-positionParticles
-{
-// ordered positioning
-	method positionOrdered;     
-// maximum number of particles in the simulation
-	maxNumberOfParticles 40000;
-// perform initial sorting based on morton code? 
-	mortonSorting Yes;             
-// cylinder for positioning particles 
-	cylinder
-	{
-// Coordinates of top cylinderRegion (m,m,m)	
-		p1 (0.05 0.0 0.12);
-		p2 (0.05 0.0 0.22);
-// radius of cylinder
-		radius 0.066;
-	}
-
-	positionOrderedInfo
-	{
-// minimum space between centers of particles
-		diameter 0.003;
-// number of particles in the simulation 	 	
-		numPoints 20000;
-// axis order for filling the space with particles		 	
-		axisOrder (z y x);  
-	}
-}
-```
- ## Interaction between particles
- In `caseSetup/interaction` file, material names and properties and interaction parameters are defined: interaction between the particles of rotating drum. Since we are defining 1 material for simulation, the interaction matrix is 1x1 (interactions are symetric). 
-```C++
- // a list of materials names
-materials      (prop1);
-// density of materials [kg/m3]
-densities      (1000.0);   
-
-contactListType   sortedContactList; 
-
-model
-{
-   contactForceModel nonLinearNonLimited;
-   rollingFrictionModel normal;
-   /*
-   Property (prop1-prop1);
-   */
-// Young modulus [Pa]
-   Yeff  (1.0e6);       
-// Shear modulus [Pa]
-   Geff  (0.8e6);       
-// Poisson's ratio [-]
-   nu    (0.25);        
-// coefficient of normal restitution
-   en    (0.7);         
-// coefficient of tangential restitution
-   et    (1.0);          
-// dynamic friction
-   mu    (0.3);          
-// rolling friction
-   mur   (0.1);          
-        
-}
-```
-## Settings
-### Geometry
-In the `settings/geometryDict` file, the geometry and axis of rotation is defined for the drum. The geometry is composed of a cylinder inlet and outlet, cone shell top and down, a cylinder shell and enter and exit Gate.
-```C++
-surfaces
-{
-	topGate
-	topGate
-	{
-	// type of wall
-		type cylinderWall;
-	// begin point of cylinder axis 
-		p1 (0.0    0.0   0.299);
-	// end point of cylinder axis 
-		p2 (0.0    0.0   0.3);
-	// radius at p1  
-		radius1  0.03;
-	// radius at p2		
-		radius2	 0.0001;
-	// material of wall
-		material solidProperty;
-	// motion component name
-		motion axisOfRotation;	
-	}
-
-    topCylinder
-	{
-	// type of the wall
-		type 		cylinderWall;
-	// begin point of cylinder axis 	
-		p1 			(0.0 0.0 0.28);
-	// end point of cylinder axis  
-		p2 			(0.0 0.0 0.3);
-	// radius at p1  
-		radius1 	0.03;
-	// radius at p2		
-		radius2 	0.03;
-	// number of divisions		
-		resolution 	36;
-	// material name of this wall
-		material 	prop1;
-	// motion component name   	
-		motion axisOfRotation;		
-	}
-
-	coneShelltop
-	{
-	// type of the wall	
-		type 		cylinderWall;
-	// begin point of cylinder axis  	
-		p1 			(0.0 0.0 0.2);
-	// end point of cylinder axis  
-		p2 			(0.0 0.0 0.28);
-	// radius at p1  
-		radius1 	0.1;
-	// radius at p2		
-		radius2 	0.03;
-	// number of divisions		
-		resolution 	36;
-	// material name of this wall      	
-		material 	prop1;
-	// motion component name   	
-		motion axisOfRotation;		
-	}
-
-	cylinderShell
-	{
-	// type of the wall
-		type 		cylinderWall;
-	// begin point of cylinder axis  	
-		p1 			(0.0 0.0 0.1);
-	// end point of cylinder axis
-		p2 			(0.0 0.0 0.2);
-	// radius at p1	  
-		radius1 	0.1;
-	// radius at p2			
-		radius2 	0.1;
-	// number of divisions		
-		resolution 	36;
-	// material name of this wall	      	
-		material 	prop1; 
-	// motion component name  	
-		motion axisOfRotation;		
-	}
-
-	coneShelldown
-	{
-	// type of the wall
-		type 		cylinderWall;
-	// begin point of cylinder axis  	
-		p1 			(0.0 0.0 0.02);
-	// end point of cylinder axis  
-		p2 			(0.0 0.0 0.1);
-	// radius at p1  
-		radius1 	0.03;
-	// radius at p2		
-		radius2 	0.1;
-	// number of divisions		
-		resolution 	36;
-	// material name of this wall	      	
-		material 	prop1;
-	// motion component name   	
-		motion axisOfRotation;		
-	}
-	/*
-	This is a plane wall at the exit of silo
-	*/
-
-	    bottomCylinder
-	{
-	// type of the wall
-		type 		cylinderWall;  	
-	// begin point of cylinder axis	
-		p1 			(0.0 0.0 0.0);
-	// end point of cylinder axis	  
-		p2 			(0.0 0.0 0.02);
-	// radius at p1  
-		radius1 	0.03;
-	// radius at p2			
-		radius2 	0.03;
-	// number of divisions			
-		resolution 	36;
-	// material name of this wall	      	
-		material 	prop1;
-	// motion component name	   	
-		motion axisOfRotation;		
-	}
-	exitGate
-	{
-		type planeWall;
-		p1 (-0.05    -0.05    0);
-		p2 (-0.05    0.05     0);
-		p3 ( 0.05    0.05     0);
-		p4 (0.05     -0.05    0);
-		material prop1;
-		motion axisOfRotation;		
-	}
-		
-}
-```
-### Rotating Axis Info
-In this part of `geometryDict` the information of rotating axis and speed of rotation are defined. Unlike the previous cases, the rotation of this blender starts at time=**0 s**.
 ```C++
 // information for rotatingAxisMotion motion model 
 rotatingAxisMotionInfo
@@ -248,18 +32,325 @@ rotatingAxisMotionInfo
 	{
 		p1 (-0.1 0.0 0.15);	// first point for the axis of rotation 
 		p2 ( 0.1 0.0 0.15);	// second point for the axis of rotation
+		
 		omega 1.5708; 		// rotation speed ==> 15 rad/s
+		
 		// Start time of Geometry Rotating (s) 		
-		startTime 1;
+		startTime 0.5;
+		
 		// End time of Geometry Rotating (s)
 		endTime 9.5;
 	}
 }
 ```
-## Performing Simulation
+
+
+### Surfaces 
+In `settings/geometryDict` file, the surfaces and motion component of each surface are defined to form a rotating tote-blender. The geometry is composed of top and bottom cylinders, top and bottom cones, a cylindrical shell and top and bottom Gates.
+
+```C++
+surfaces
+{
+	
+	topGate
+	{
+		// type of wall
+		type cylinderWall;
+		
+		// begin point of cylinder axis 
+		p1 (0.0    0.0   0.3);
+		
+		// end point of cylinder axis 
+		p2 (0.0    0.0   0.301);
+		
+		// radius at p1  
+		radius1  0.03;
+		
+		// radius at p2		
+		radius2	 0.0001;
+		
+		// material of wall
+		material solidProperty;
+		
+		// motion component name
+		motion axisOfRotation;	
+	}
+
+    topCylinder
+	{
+		// type of the wall
+		type 		cylinderWall;
+		
+		// begin point of cylinder axis 	
+		p1 			(0.0 0.0 0.28);
+		
+		// end point of cylinder axis  
+		p2 			(0.0 0.0 0.3);
+		
+		// radius at p1  
+		radius1 	0.03;
+		
+		// radius at p2		
+		radius2 	0.03;
+		
+		// number of divisions		
+		resolution 	36;
+		
+		// material name of this wall
+		material 	solidProperty;
+		
+		// motion component name   	
+		motion axisOfRotation;		
+	}
+
+	coneShelltop
+	{	
+		// type of the wall	
+		type 		cylinderWall;
+		
+		// begin point of cylinder axis  	
+		p1 			(0.0 0.0 0.2);
+		
+		// end point of cylinder axis  
+		p2 			(0.0 0.0 0.28);
+		
+		// radius at p1  
+		radius1 	0.1;
+		
+		// radius at p2		
+		radius2 	0.03;
+		
+		// number of divisions		
+		resolution 	36;
+		
+		// material name of this wall      	
+		material 	solidProperty;
+		
+		// motion component name   	
+		motion axisOfRotation;		
+	}
+
+	cylinderShell
+	{
+		// type of the wall
+		type 		cylinderWall;
+		
+		// begin point of cylinder axis  	
+		p1 			(0.0 0.0 0.1);
+		
+		// end point of cylinder axis
+		p2 			(0.0 0.0 0.2);
+		
+		// radius at p1	  
+		radius1 	0.1;
+		
+		// radius at p2			
+		radius2 	0.1;
+		
+		// number of divisions		
+		resolution 	36;
+		
+		// material name of this wall	      	
+		material 	solidProperty; 
+		
+		// motion component name  	
+		motion axisOfRotation;		
+	}
+
+	coneShelldown
+	{
+		
+		// type of the wall
+		type 		cylinderWall;
+		
+		// begin point of cylinder axis  	
+		p1 			(0.0 0.0 0.02);
+		
+		// end point of cylinder axis  
+		p2 			(0.0 0.0 0.1);
+		
+		// radius at p1  
+		radius1 	0.03;
+		
+		// radius at p2		
+		radius2 	0.1;
+		
+		// number of divisions		
+		resolution 	36;
+		
+		// material name of this wall	      	
+		material 	solidProperty;
+		
+		// motion component name   	
+		motion axisOfRotation;		
+	}
+
+	bottomCylinder
+	{
+		// type of the wall
+		type 		cylinderWall;  	
+		
+		// begin point of cylinder axis	
+		p1 			(0.0 0.0 0.0);
+		
+		// end point of cylinder axis	  
+		p2 			(0.0 0.0 0.02);
+		
+		// radius at p1  
+		radius1 	0.03;
+		
+		// radius at p2			
+		radius2 	0.03;
+		
+		// number of divisions			
+		resolution 	36;
+		
+		// material name of this wall	      	
+		material 	solidProperty;
+		
+		// motion component name	   	
+		motion axisOfRotation;		
+	}
+
+	exitGate
+	{
+		
+		// type of the wall
+		type 		cylinderWall;  	
+		
+		// begin point of cylinder axis	
+		p1 			(0.0 0.0 -0.001);
+		
+		// end point of cylinder axis	  
+		p2 			(0.0 0.0 0.0);
+		
+		// radius at p1  
+		radius1 	0.03;
+		
+		// radius at p2			
+		radius2 	0.0001;
+		
+		// number of divisions			
+		resolution 	36;
+		
+		// material name of this wall	      	
+		material 	solidProperty;
+		
+		// motion component name	   	
+		motion axisOfRotation;			
+	}
+		
+}
+```
+
+## Defining particles 
+### Diameter and material of spheres 
+In the `caseSetup/sphereShape` the diameter and the material name of the particles are defined.
+
+<div align="center"> 
+in <b>caseSetup/sphereShape</b> file
+</div>
+
+```C++
+// name of shapes 
+names 		(sphere1); 	
+
+// diameter of shapes (m)
+diameters 	(0.005);
+
+// material name for shapes 	
+materials	(solidProperty);
+```
+### Particle positioning before start of simulation 
+Particles are positioned before the start of simulation. The positioning can be ordered or random. Here we use ordered positioning. 24000 particles are positioned in a cylinderical region inside the tote-blender.
+
+<div align="center"> 
+in <b>settings/particlesDict</b> file
+</div>
+
+```C++
+// positions particles 
+positionParticles
+{
+	// ordered positioning
+	method 	positionOrdered;     
+
+	// maximum number of particles in the simulation
+	maxNumberOfParticles 25001;
+	
+	// perform initial sorting based on morton code? 
+	mortonSorting 	Yes;             
+	
+	// cylinderical region for positioning particles 
+	cylinder
+	{
+		p1 (0.0 0.0 0.09);	
+		p2 (0.0 0.0 0.21);
+		radius 0.09;
+	}
+
+	positionOrderedInfo
+	{
+		// minimum space between centers of particles
+		diameter 0.005;
+		
+		// number of particles in the simulation 	 	
+		numPoints 24000;
+	
+		// axis order for filling the space with particles		 	
+		axisOrder (x y z);  
+	}
+}
+```
+
+ ## Interaction between particles
+ In `caseSetup/interaction` file, material names and properties and interaction parameters are defined. Since we are defining 1 material type in the simulation, the interaction matrix is 1x1 (interactions are symmetric). 
+```C++
+ // a list of materials names
+materials      (solidProperty);
+
+// density of materials [kg/m3]
+densities      (1000.0);   
+
+contactListType   sortedContactList; 
+
+model
+{
+   contactForceModel       nonLinearNonLimited;
+    
+   rollingFrictionModel    normal;
+   
+   /*
+   Property (solidProperty-solidProperty);
+   */
+    
+   // Young modulus [Pa]
+   Yeff  (1.0e6);       
+
+   // Shear modulus [Pa]
+   Geff  (0.8e6);       
+    
+   // Poisson's ratio [-]
+   nu    (0.25);        
+    
+   // coefficient of normal restitution
+   en    (0.7);         
+
+   // coefficient of tangential restitution
+   et    (1.0);          
+
+   // dynamic friction
+   mu    (0.3);          
+
+   // rolling friction
+   mur   (0.1);                  
+}
+```
+
+# Performing Simulation and previewing the results 
 To perform simulations, enter the following commands one after another in the terminal. 
 
 Enter `$ particlesPhasicFlow` command to create the initial fields for particles.  
-Enter `$ geometryPhasicFlow` command to create the Geometry.  
+Enter `$ geometryPhasicFlow` command to create the geometry.  
 At last, enter `$ sphereGranFlow` command to start the simulation.  
-After finishing the simulation, you can use  `$ pFlowtoVTK` to convert the results into vtk format storred in ./VTK folder.
+After finishing the simulation, you can use  `$ pFlowtoVTK` to convert the results into vtk format stored in ./VTK folder.