From 7f3b4e0e2cff635b67c7ee49c64d7b3ee6f31111 Mon Sep 17 00:00:00 2001
From: hamidrezanorouzi <hamid.r.norouzi@gmail.com>
Date: Mon, 5 Sep 2022 11:14:41 +0430
Subject: [PATCH] utilities added

---
 .gitignore                                    |   1 +
 CMakeLists.txt                                |   2 +-
 utilities/CMakeLists.txt                      |  10 +
 utilities/geometryPhasicFlow/CMakeLists.txt   |  12 +
 utilities/geometryPhasicFlow/Wall/Wall.C      | 103 +++++
 utilities/geometryPhasicFlow/Wall/Wall.H      | 113 ++++++
 .../cylinderWall/cylinderWall.C               | 101 +++++
 .../cylinderWall/cylinderWall.H               |  63 +++
 .../geometryPhasicFlow/cylinderWall/zAxis.C   | 150 ++++++++
 .../geometryPhasicFlow/cylinderWall/zAxis.H   | 107 ++++++
 .../geometryPhasicFlow/geometryPhasicFlow.C   |  98 +++++
 .../geometryPhasicFlow/planeWall/planeWall.C  |  75 ++++
 .../geometryPhasicFlow/planeWall/planeWall.H  |  60 +++
 .../geometryPhasicFlow/stlWall/stlWall.C      |  71 ++++
 .../geometryPhasicFlow/stlWall/stlWall.H      |  60 +++
 utilities/pFlowToVTK/CMakeLists.txt           |   9 +
 utilities/pFlowToVTK/geometric.C              |  87 +++++
 utilities/pFlowToVTK/geometric.H              |  78 ++++
 utilities/pFlowToVTK/pFlowToVTK.C             | 142 +++++++
 utilities/pFlowToVTK/pointFieldToVTK.H        | 361 ++++++++++++++++++
 utilities/pFlowToVTK/triSurfaceFieldToVTK.H   | 245 ++++++++++++
 utilities/pFlowToVTK/vtkFile.C                |  75 ++++
 utilities/pFlowToVTK/vtkFile.H                |  90 +++++
 utilities/particlesPhasicFlow/CMakeLists.txt  |  11 +
 utilities/particlesPhasicFlow/empty/empty.C   |  36 ++
 utilities/particlesPhasicFlow/empty/empty.H   |  94 +++++
 .../particlesPhasicFlow/particlesPhasicFlow.C | 195 ++++++++++
 .../positionOrdered/positionOrdered.C         | 161 ++++++++
 .../positionOrdered/positionOrdered.H         | 115 ++++++
 .../positionParticles/positionParticles.C     | 110 ++++++
 .../positionParticles/positionParticles.H     |  98 +++++
 .../positionRandom/positionRandom.C           | 260 +++++++++++++
 .../positionRandom/positionRandom.H           | 113 ++++++
 utilities/particlesPhasicFlow/setFields.H     |  56 +++
 34 files changed, 3361 insertions(+), 1 deletion(-)
 create mode 100644 utilities/CMakeLists.txt
 create mode 100644 utilities/geometryPhasicFlow/CMakeLists.txt
 create mode 100755 utilities/geometryPhasicFlow/Wall/Wall.C
 create mode 100755 utilities/geometryPhasicFlow/Wall/Wall.H
 create mode 100644 utilities/geometryPhasicFlow/cylinderWall/cylinderWall.C
 create mode 100644 utilities/geometryPhasicFlow/cylinderWall/cylinderWall.H
 create mode 100644 utilities/geometryPhasicFlow/cylinderWall/zAxis.C
 create mode 100644 utilities/geometryPhasicFlow/cylinderWall/zAxis.H
 create mode 100755 utilities/geometryPhasicFlow/geometryPhasicFlow.C
 create mode 100755 utilities/geometryPhasicFlow/planeWall/planeWall.C
 create mode 100755 utilities/geometryPhasicFlow/planeWall/planeWall.H
 create mode 100755 utilities/geometryPhasicFlow/stlWall/stlWall.C
 create mode 100755 utilities/geometryPhasicFlow/stlWall/stlWall.H
 create mode 100644 utilities/pFlowToVTK/CMakeLists.txt
 create mode 100755 utilities/pFlowToVTK/geometric.C
 create mode 100755 utilities/pFlowToVTK/geometric.H
 create mode 100755 utilities/pFlowToVTK/pFlowToVTK.C
 create mode 100755 utilities/pFlowToVTK/pointFieldToVTK.H
 create mode 100644 utilities/pFlowToVTK/triSurfaceFieldToVTK.H
 create mode 100755 utilities/pFlowToVTK/vtkFile.C
 create mode 100755 utilities/pFlowToVTK/vtkFile.H
 create mode 100644 utilities/particlesPhasicFlow/CMakeLists.txt
 create mode 100755 utilities/particlesPhasicFlow/empty/empty.C
 create mode 100755 utilities/particlesPhasicFlow/empty/empty.H
 create mode 100755 utilities/particlesPhasicFlow/particlesPhasicFlow.C
 create mode 100755 utilities/particlesPhasicFlow/positionOrdered/positionOrdered.C
 create mode 100755 utilities/particlesPhasicFlow/positionOrdered/positionOrdered.H
 create mode 100755 utilities/particlesPhasicFlow/positionParticles/positionParticles.C
 create mode 100755 utilities/particlesPhasicFlow/positionParticles/positionParticles.H
 create mode 100755 utilities/particlesPhasicFlow/positionRandom/positionRandom.C
 create mode 100755 utilities/particlesPhasicFlow/positionRandom/positionRandom.H
 create mode 100755 utilities/particlesPhasicFlow/setFields.H

diff --git a/.gitignore b/.gitignore
index ed61ef2a..d0570535 100644
--- a/.gitignore
+++ b/.gitignore
@@ -36,3 +36,4 @@ build/**
 include/**
 bin/**
 lib/**
+**/**notnow
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e34b2463..3b3651a9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -72,7 +72,7 @@ add_subdirectory(src)
 
 add_subdirectory(solvers)
 
-#add_subdirectory(utilities)
+add_subdirectory(utilities)
 
 #add_subdirectory(test)
 
diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
new file mode 100644
index 00000000..1be15174
--- /dev/null
+++ b/utilities/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+
+
+add_subdirectory(particlesPhasicFlow)
+
+add_subdirectory(geometryPhasicFlow)
+
+add_subdirectory(pFlowToVTK)
+
+
diff --git a/utilities/geometryPhasicFlow/CMakeLists.txt b/utilities/geometryPhasicFlow/CMakeLists.txt
new file mode 100644
index 00000000..31aed235
--- /dev/null
+++ b/utilities/geometryPhasicFlow/CMakeLists.txt
@@ -0,0 +1,12 @@
+
+set(source_files 
+geometryPhasicFlow.C
+Wall/Wall.C 
+planeWall/planeWall.C 
+stlWall/stlWall.C
+cylinderWall/zAxis.C
+cylinderWall/cylinderWall.C
+)
+set(link_lib phasicFlow Geometry Kokkos::kokkos)
+
+pFlow_make_executable_install(geometryPhasicFlow source_files link_lib)
diff --git a/utilities/geometryPhasicFlow/Wall/Wall.C b/utilities/geometryPhasicFlow/Wall/Wall.C
new file mode 100755
index 00000000..9c5f539e
--- /dev/null
+++ b/utilities/geometryPhasicFlow/Wall/Wall.C
@@ -0,0 +1,103 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#include "Wall.H"
+
+bool pFlow::Wall::readCommon(const dictionary& dict)
+{
+	materialName_ = dict.getVal<word>("material");
+
+	motionName_ = dict.getValOrSet("motion", word("none"));
+
+	name_ = dict.name();
+
+	return true;
+}
+
+
+pFlow::Wall::Wall(const dictionary& dict)
+{
+	if(!readCommon(dict))
+	{
+		fatalExit;
+	}
+}
+
+
+bool pFlow::Wall::checkTrianlge
+(
+	const realx3 &p1,
+	const realx3 &p2,
+	const realx3 &p3
+)
+{
+	realx3 ln = cross(p2 - p1, p3 - p1);
+	
+	if (ln.length() < smallValue) return false;
+	
+	return true;
+}
+
+pFlow::uniquePtr<pFlow::Wall> 
+pFlow::Wall::create
+(
+	const dictionary& dict
+)
+{
+	word type = dict.getVal<word>("type");
+
+	if( dictionaryvCtorSelector_.search(type) )
+	{
+		return dictionaryvCtorSelector_[type] (dict);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< type << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+	return nullptr;
+}
+
+namespace pFlow
+{
+
+bool checkNormalVec(
+	const realx3 &p1,
+	const realx3 &p2,
+	const realx3 &p3,
+	realx3& norm )
+{
+	realx3 ln = cross(p2 - p1, p3 - p1);
+	real len = length(ln);
+	if (len < smallValue) return false;
+	
+	norm = ln/len;
+
+	return true;
+}
+
+}
\ No newline at end of file
diff --git a/utilities/geometryPhasicFlow/Wall/Wall.H b/utilities/geometryPhasicFlow/Wall/Wall.H
new file mode 100755
index 00000000..d08a1c6e
--- /dev/null
+++ b/utilities/geometryPhasicFlow/Wall/Wall.H
@@ -0,0 +1,113 @@
+/*------------------------------- 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 __Wall_H__
+#define __Wall_H__
+
+
+#include "virtualConstructor.H"
+#include "Vectors.H"
+#include "dictionary.H"
+
+namespace pFlow
+{
+
+
+bool checkNormalVec(
+	const realx3 &p1,
+	const realx3 &p2,
+	const realx3 &p3,
+	realx3& norm );
+
+
+class Wall
+{
+protected:
+
+	realx3x3Vector 	triangles_;
+
+	word   			name_;
+
+	word  			materialName_;
+
+	word  			motionName_;
+
+	bool readCommon(const dictionary& ditc);
+
+public:
+
+	// - type info
+	TypeName("Wall");
+
+	//// - Constructors
+
+		// - empty
+		Wall(){}
+
+		// - from dictionary
+		Wall(const dictionary& dict);
+
+		
+		virtual ~Wall()=default;
+
+		create_vCtor
+		(
+			Wall,
+			dictionary,
+			(const dictionary& dict),
+			(dict)
+		);
+
+	//// - Methods
+
+		// - 
+		const realx3x3Vector& triangles()const
+		{
+			return triangles_;
+		}
+
+		word name()const
+		{
+			return name_;
+		}
+
+		word materialName()const
+		{
+			return materialName_;
+		}
+
+		word motionName()const
+		{
+			return motionName_;
+		}
+
+		static 
+			bool checkTrianlge(const realx3& p1, const realx3& p2, const realx3& p3);
+
+		static
+			uniquePtr<Wall> create(const dictionary& dict);
+
+};
+
+}
+
+
+#endif
diff --git a/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.C b/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.C
new file mode 100644
index 00000000..bf10a5d8
--- /dev/null
+++ b/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.C
@@ -0,0 +1,101 @@
+#include "cylinderWall.H"
+
+
+bool pFlow::cylinderWall::readCylinderWall(const dictionary& dict)
+{
+	auto p1 = dict.getVal<realx3>("p1");
+	auto p2 = dict.getVal<realx3>("p2");
+	auto radius1 = dict.getVal<real>("radius1");
+	auto radius2 = dict.getVal<real>("radius2") ;
+
+	int32 resolution = dict.getValOrSet("resolution", 24 );
+
+
+	return createCylinder(p1, p2, radius1, radius2, resolution);
+}
+
+bool pFlow::cylinderWall::createCylinder(
+	const realx3& p1,
+	const realx3& p2,
+	real rad1,
+	real rad2,
+	int32 numDiv)
+{
+
+	zAxis zAx(p1, p2);
+	
+	real L = zAx.length();
+
+	// number of wall elements will be twice numDiv
+	triangles_.clear();
+	triangles_.reserve(2 * numDiv);
+
+	realx3Vector r1P(numDiv + 1), r2P(numDiv + 1);
+	real dTheta = 2 * Pi / numDiv;
+	real theta = 0;
+
+	for (int32 i = 0; i < numDiv + 1; i++)
+	{
+		r1P[i] = realx3(rad1*cos(theta), rad1*sin(theta), 0);
+		r2P[i] = realx3(rad2*cos(theta), rad2*sin(theta), L);
+		theta += dTheta;
+	}
+	
+	// transferring back all points to the original axis of cylinder
+	for (int32 i = 0; i < numDiv + 1; i++)
+	{
+		r1P[i] = zAx.transferBackZ(r1P[i]);
+		r2P[i] = zAx.transferBackZ(r2P[i]);
+	}
+	
+	realx3 norm;
+	for (int32 i = 0; i < numDiv; i++)
+	{
+		realx3 p1 = r1P[i];
+		realx3 p2 = r2P[i];
+		realx3 p3 = r2P[i + 1];
+		realx3 p4 = r1P[i + 1];
+
+		if(checkNormalVec(p1, p2, p3, norm))
+		{
+			triangles_.push_back(realx3x3(p1, p2, p3));
+		}
+		else
+		{
+			fatalErrorInFunction<<
+			"planner input triangle: "<<realx3x3(p1, p2, p3)<<endl;
+			return false;
+		}
+		
+		if (checkNormalVec(p3, p4, p1, norm))
+		{
+			triangles_.push_back(realx3x3(p3, p4, p1));
+		}
+		else
+		{
+			fatalErrorInFunction<<
+			"planner input triangle: "<<realx3x3(p3, p4, p1)<<endl;
+			return false;
+		}
+
+	}
+
+	return true;
+}
+
+
+pFlow::cylinderWall::cylinderWall()
+{}
+
+
+pFlow::cylinderWall::cylinderWall(
+	const dictionary& dict
+)
+:
+	Wall(dict)
+{
+	if( !readCylinderWall(dict) )
+	{
+		fatalExit;
+	}
+}
\ No newline at end of file
diff --git a/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.H b/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.H
new file mode 100644
index 00000000..e8533aa4
--- /dev/null
+++ b/utilities/geometryPhasicFlow/cylinderWall/cylinderWall.H
@@ -0,0 +1,63 @@
+/*------------------------------- 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 __cylinderWall_H__
+#define __cylinderWall_H__
+
+#include "Wall.H"
+#include "zAxis.H"
+
+namespace pFlow
+{
+
+class cylinderWall
+:
+	public Wall
+{
+protected:
+
+	bool readCylinderWall(const dictionary& dict);
+
+	bool createCylinder(const realx3& p1, const realx3& p2, real rad1, real rad2, int32 numDiv);
+
+public:
+
+	TypeName("cylinderWall");
+
+	cylinderWall();
+
+	cylinderWall(const dictionary& dict);
+
+	add_vCtor
+	(
+		Wall,
+		cylinderWall,
+		dictionary
+	);
+
+};
+
+} // pFlow
+
+
+
+
+#endif //__cylinderWall_H__
\ No newline at end of file
diff --git a/utilities/geometryPhasicFlow/cylinderWall/zAxis.C b/utilities/geometryPhasicFlow/cylinderWall/zAxis.C
new file mode 100644
index 00000000..af1b0725
--- /dev/null
+++ b/utilities/geometryPhasicFlow/cylinderWall/zAxis.C
@@ -0,0 +1,150 @@
+/*------------------------------- 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 "zAxis.H"
+
+
+pFlow::zAxis::zAxis(const realx3 &p1, const realx3 &p2) :
+p1_(p1),
+p2_(p2)
+{
+	n_ = p2-p1;
+	auto len = pFlow::length(n_);
+
+	if(len < smallValue )
+	{
+		fatalExit;
+	}
+	n_ /= len;
+
+	makeTransMatrix();
+}
+
+pFlow::realx3 pFlow::zAxis::transferToZ(const realx3 & p)
+{
+	real pp[4][1] = {p.x(), p.y(), p.z(), 1.0};
+	real pn[4][1];
+
+	MatMul(Trans_z_xz_P1_, pp, pn);
+	
+	return realx3(
+		pn[0][0],
+		pn[1][0],
+		pn[2][0]
+		);
+
+}
+
+pFlow::realx3 pFlow::zAxis::transferBackZ(const realx3 & p)
+{
+	real pp[4][1] = { p.x(), p.y(), p.z(), 1.0 };
+	real pn[4][1];
+
+	MatMul(ITrans_P1_xz_z_, pp, pn);
+
+	return realx3(
+		pn[0][0],
+		pn[1][0],
+		pn[2][0]
+		);
+}
+
+
+void pFlow::zAxis::makeTransMatrix()
+{
+	
+	// transfering point p1 to the origin
+	real TransP1[4][4] =
+	{
+		1.0, 0.0, 0.0, -p1_.x(),
+		0.0, 1.0, 0.0, -p1_.y(),
+		0.0, 0.0, 1.0, -p1_.z(),
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	// for transformation back
+	real ITransP1[4][4] =
+	{
+		1.0, 0.0, 0.0, p1_.x(),
+		0.0, 1.0, 0.0, p1_.y(),
+		0.0, 0.0, 1.0, p1_.z(),
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	
+	real u = n_.x();
+	real v = n_.y();
+	real w = n_.z();
+	real u2v2 = sqrt(u*u + v*v);
+
+	//correcting the transformation matrix in the case of coincidence with z - axis
+	if ( equal(w,1.0) )
+	{
+		assignMat(TransP1 , Trans_z_xz_P1_);
+		assignMat(ITransP1, ITrans_P1_xz_z_);
+		return;
+	}
+
+	u2v2 = max(smallValue, u2v2);
+
+	real TransXZ[4][4] =
+	{
+		u / u2v2, v / u2v2, 0.0, 0.0,
+		-v / u2v2, u / u2v2, 0.0, 0.0,
+		0.0, 0.0, 1.0, 0.0,
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	real TransZ[4][4] =
+	{
+		w, 0.0, -u2v2, 0.0,
+		0.0, 1.0, 0.0, 0.0,
+		u2v2, 0.0, w, 0.0,
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	real temp[4][4];
+	// creat transformation matrix to transfer point from line axis to z-axis
+	MatMul(TransXZ, TransP1, temp);
+	MatMul(TransZ, temp, Trans_z_xz_P1_);
+
+
+	real ITransXZ[4][4] =
+	{
+		u / u2v2, -v / u2v2, 0.0, 0.0,
+		+v / u2v2, u / u2v2, 0.0, 0.0,
+		0.0, 0.0, 1.0, 0.0,
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	real ITransZ[4][4] =
+	{
+		w, 0.0, +u2v2, 0.0,
+		0.0, 1.0, 0.0, 0.0,
+		-u2v2, 0.0, w, 0.0,
+		0.0, 0.0, 0.0, 1.0
+	};
+
+	// creat transformation matrix to transfer point to from z-axis to line axis
+	MatMul(ITransXZ, ITransZ, temp);
+	MatMul(ITransP1, temp, ITrans_P1_xz_z_);
+
+}
diff --git a/utilities/geometryPhasicFlow/cylinderWall/zAxis.H b/utilities/geometryPhasicFlow/cylinderWall/zAxis.H
new file mode 100644
index 00000000..d0a455b0
--- /dev/null
+++ b/utilities/geometryPhasicFlow/cylinderWall/zAxis.H
@@ -0,0 +1,107 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+/*
+class name: zAxis
+a class for transfering points to a cylinderical coordinates defined by points
+p1 and p2
+*/
+
+#ifndef __zAxis_H__
+#define __zAxis_H__
+
+#include "types.H"
+
+namespace pFlow
+{
+
+template <typename T, int32 nRow, int32 nInner, int32 nCol >
+void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol]);
+
+template <typename T, int32 nRow, int32 nCol >
+void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol]);
+
+class zAxis
+{
+public:
+	// constructors 
+	zAxis(const realx3 &lp1, const realx3 &lp2);
+	
+	real length()const
+	{
+		return pFlow::length(p2_-p1_);
+	}
+
+	realx3 transferToZ(const realx3 & p);
+
+	realx3 transferBackZ(const realx3 & p);
+
+private:
+	void makeTransMatrix();
+protected:
+
+	realx3 p1_;
+	realx3 p2_;
+	realx3 n_;
+
+	real Trans_z_xz_P1_[4][4];
+
+	real ITrans_P1_xz_z_[4][4];
+};
+
+
+
+template <typename T, int32 nRow, int32 nInner, int32 nCol >
+void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol])
+{
+
+	for (int32 row = 0; row < nRow; row++)
+	{
+		for (int32 col = 0; col < nCol; col++)
+		{
+			T sum = 0;
+			for (int inner = 0; inner < nInner; inner++)
+			{
+				sum += A[row][inner] * B[inner][col];
+			}
+			C[row][col] = sum;
+		}
+	}
+}
+
+template <typename T, int32 nRow, int32 nCol >
+void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol])
+{
+
+	for (int32 row = 0; row < nRow; row++)
+	{
+		for (int32 col = 0; col < nCol; col++)
+		{
+			B[row][col] = A[row][col];
+		}
+	}
+}
+
+
+}
+
+#endif
+
diff --git a/utilities/geometryPhasicFlow/geometryPhasicFlow.C b/utilities/geometryPhasicFlow/geometryPhasicFlow.C
new file mode 100755
index 00000000..26f03bf0
--- /dev/null
+++ b/utilities/geometryPhasicFlow/geometryPhasicFlow.C
@@ -0,0 +1,98 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#include "systemControl.H"
+#include "Wall.H"
+#include "multiTriSurface.H"
+#include "geometryMotion.H"
+
+using pFlow::output;
+using pFlow::endl;
+using pFlow::IOobject;
+using pFlow::dictionary;
+using pFlow::objectFile;
+using pFlow::wordVector;
+using pFlow::Wall;
+using pFlow::geometry;
+using pFlow::multiTriSurface;
+
+int main( int argc, char* argv[] )
+{
+
+// this should be palced in each main 
+#include "initialize_Control.H"
+
+	#include "setProperty.H"
+
+	Report(0)<<"\nReading "<<"createGeometryDict"<<" . . ."<<endReport;
+	auto objDict = IOobject::make<dictionary>
+	(
+		objectFile
+		(
+			"createGeometryDict",
+			Control.settings().path(),
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_NEVER
+		),
+		"createGeometryDict",
+		true
+	);
+
+	auto& geometryDict = objDict().getObject<dictionary>();
+
+	auto& surfacesDict = geometryDict.subDict("surfaces");
+
+	auto wallsDictName = surfacesDict.dictionaryKeywords();
+
+	
+
+	multiTriSurface surface;
+	wordVector materials;
+	wordVector motion;
+	
+	for(auto& name:wallsDictName)
+	{
+		Report(1)<<"Creating wall "<<greenText(name)<<" from file dictionary . . . "<<endReport;
+		auto wallPtr = Wall::create( surfacesDict.subDict(name));
+		auto& wall = wallPtr();
+		Report(1)<<"wall type is "<<greenText(wall.typeName())<<'\n'<<endReport;		
+
+		surface.addTriSurface(wall.name(), wall.triangles());
+		materials.push_back(wall.materialName());
+		motion.push_back(wall.motionName());
+	}
+
+	Report(1)<<"Selected wall materials are "<<cyanText(materials)<<'\n'<<endReport;
+		
+	Report(0)<< "\nCreating geometry . . ."<<endReport;
+	auto geomPtr = geometry::create(Control, proprties, geometryDict, surface, motion, materials);
+	Report(1)<< "geometry type is "<< greenText(geomPtr().typeName())<<endReport;
+	
+	Report(1)<< "Writing geometry to folder "<< geomPtr().path()<<endl;
+	geomPtr().write();
+
+	Report(0)<< greenText("\nFinished successfully.\n");
+
+// this should be palced in each main 
+#include "finalize.H"
+
+	return 0;
+} 
diff --git a/utilities/geometryPhasicFlow/planeWall/planeWall.C b/utilities/geometryPhasicFlow/planeWall/planeWall.C
new file mode 100755
index 00000000..64064b8f
--- /dev/null
+++ b/utilities/geometryPhasicFlow/planeWall/planeWall.C
@@ -0,0 +1,75 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#include "planeWall.H"
+
+
+
+bool pFlow::planeWall::readPlaneWall
+(
+	const dictionary& dict
+)
+{
+	auto p1 = dict.getVal<realx3>("p1");
+	auto p2 = dict.getVal<realx3>("p2");
+	auto p3 = dict.getVal<realx3>("p3");
+	auto p4 = dict.getVal<realx3>("p4");
+
+	
+	if( Wall::checkTrianlge(p1,p2,p3) )
+	{
+		triangles_.push_back(realx3x3(p1,p2,p3));
+	}else
+	{
+		fatalErrorInFunction << 
+		"points p1, p2 and p3 do not form a plane wall in dictionary " << dict.globalName()<<endl;
+		return false;
+	}
+
+	if( Wall::checkTrianlge(p3,p4,p1) )
+	{
+		triangles_.push_back(realx3x3(p3,p4,p1));
+	}else
+	{
+		fatalErrorInFunction << 
+		"points p3, p4 and p1 do not form a plane wall in dictionary " << dict.globalName()<<endl;
+		return false;
+	}
+
+	return true;
+	
+}
+
+pFlow::planeWall::planeWall()
+{}
+
+pFlow::planeWall::planeWall
+(
+	const dictionary& dict
+)
+:
+	Wall(dict)
+{
+	if(!readPlaneWall(dict))
+	{
+		fatalExit;
+	}
+}
diff --git a/utilities/geometryPhasicFlow/planeWall/planeWall.H b/utilities/geometryPhasicFlow/planeWall/planeWall.H
new file mode 100755
index 00000000..cfb98aed
--- /dev/null
+++ b/utilities/geometryPhasicFlow/planeWall/planeWall.H
@@ -0,0 +1,60 @@
+/*------------------------------- 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 __planeWall_H__
+#define __planeWall_H__
+
+
+#include "Wall.H"
+#include "types.H"
+
+namespace pFlow
+{
+
+class planeWall
+:
+	public Wall
+{
+protected:
+
+	bool readPlaneWall(const dictionary& dict);
+
+public:
+
+	TypeName("planeWall");
+
+	planeWall();
+
+	planeWall(const dictionary& dict);
+
+	add_vCtor
+	(
+		Wall,
+		planeWall,
+		dictionary
+	);
+
+};
+
+} // pFlow
+
+
+#endif
diff --git a/utilities/geometryPhasicFlow/stlWall/stlWall.C b/utilities/geometryPhasicFlow/stlWall/stlWall.C
new file mode 100755
index 00000000..98daae5a
--- /dev/null
+++ b/utilities/geometryPhasicFlow/stlWall/stlWall.C
@@ -0,0 +1,71 @@
+/*------------------------------- 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 "stlWall.H"
+#include "stlFile.H"
+
+#include "streams.H"
+
+bool pFlow::stlWall::readSTLWall
+(
+	const dictionary& dict
+)
+{
+	auto fileName = dict.getVal<word>("file");
+
+	
+	fileSystem file("",fileName);
+
+	stlFile stl(file);
+	if(!stl.read())
+	{
+		fatalErrorInFunction <<
+		"  error in reading stl file "<< file <<endl;
+		return false;
+	}
+
+	for(label i=0; i<stl.size(); i++)
+	{
+		auto it = triangles_.end();
+		triangles_.insert(it, stl.solid(i).begin(), stl.solid(i).end());
+	}
+
+	
+
+	return true;
+	
+}
+
+pFlow::stlWall::stlWall()
+{}
+
+pFlow::stlWall::stlWall
+(
+	const dictionary& dict
+)
+:
+	Wall(dict)
+{
+	if(!readSTLWall(dict))
+	{
+		fatalExit;
+	}
+}
diff --git a/utilities/geometryPhasicFlow/stlWall/stlWall.H b/utilities/geometryPhasicFlow/stlWall/stlWall.H
new file mode 100755
index 00000000..e14870f7
--- /dev/null
+++ b/utilities/geometryPhasicFlow/stlWall/stlWall.H
@@ -0,0 +1,60 @@
+/*------------------------------- 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 __stlWall_H__
+#define __stlWall_H__
+
+
+#include "Wall.H"
+#include "types.H"
+
+namespace pFlow
+{
+
+class stlWall
+:
+	public Wall
+{
+protected:
+
+	bool readSTLWall(const dictionary& dict);
+
+public:
+
+	TypeName("stlWall");
+
+	stlWall();
+
+	stlWall(const dictionary& dict);
+
+	add_vCtor
+	(
+		Wall,
+		stlWall,
+		dictionary
+	);
+
+};
+
+} // pFlow
+
+
+#endif
diff --git a/utilities/pFlowToVTK/CMakeLists.txt b/utilities/pFlowToVTK/CMakeLists.txt
new file mode 100644
index 00000000..5bf43fce
--- /dev/null
+++ b/utilities/pFlowToVTK/CMakeLists.txt
@@ -0,0 +1,9 @@
+
+set(source_files 
+pFlowToVTK.C 
+vtkFile.C 
+geometric.C 
+)
+set(link_lib phasicFlow Kokkos::kokkos)
+
+pFlow_make_executable_install(pFlowToVTK source_files link_lib)
diff --git a/utilities/pFlowToVTK/geometric.C b/utilities/pFlowToVTK/geometric.C
new file mode 100755
index 00000000..404e6635
--- /dev/null
+++ b/utilities/pFlowToVTK/geometric.C
@@ -0,0 +1,87 @@
+/*------------------------------- 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 "geometric.H"
+
+
+template<>
+bool pFlow::dataToVTK( vtkFile& vtk, const triSurface& surface )
+{
+	
+
+	auto nP = surface.numPoints();
+	auto hPoints = surface.points().hostVector();
+
+	vtk() << "DATASET POLYDATA" << endl;
+	vtk() << "POINTS " << nP << " float" << endl;
+
+
+	for ( auto i=0; i<nP; i++ )
+	{
+		vtk() << hPoints[i].x()  << " " << hPoints[i].y() << " " << hPoints[i].z() << endl;
+		if (!vtk) return false;
+	}
+
+	auto nV = surface.numTriangles();
+	auto hVertices = surface.vertices().hostVector();
+
+	vtk() << "POLYGONS " << nV << "  " << 4*nV << endl;
+
+	for(auto i=0; i<nV; i++)
+	{
+		vtk()<< 3 <<" "<< hVertices[i].x() << " " << hVertices[i].y() <<" "<<hVertices[i].z()<<endl;
+		if (!vtk) return false;
+	}
+
+	return true;
+} 
+
+template<>
+bool pFlow::dataToVTK( vtkFile& vtk, const multiTriSurface& surface )
+{
+
+	auto nP = surface.numPoints();
+	auto hPoints = surface.points().hostVector();
+
+	vtk() << "DATASET POLYDATA" << endl;
+	vtk() << "POINTS " << nP << " float" << endl;
+
+
+	for ( auto i=0; i<nP; i++ )
+	{
+		vtk() << hPoints[i].x()  << " " << hPoints[i].y() << " " << hPoints[i].z() << endl;
+		if (!vtk) return false;
+	}
+
+	auto nV = surface.numTriangles();
+	auto hVertices = surface.vertices().hostVector();
+
+	vtk() << "POLYGONS " << nV << "  " << 4*nV << endl;
+
+	for(auto i=0; i<nV; i++)
+	{
+		vtk()<< 3 <<" "<< hVertices[i].x() << " " << hVertices[i].y() <<" "<<hVertices[i].z()<<endl;
+		if (!vtk) return false;
+	}
+
+	return true;
+
+}
diff --git a/utilities/pFlowToVTK/geometric.H b/utilities/pFlowToVTK/geometric.H
new file mode 100755
index 00000000..13330c38
--- /dev/null
+++ b/utilities/pFlowToVTK/geometric.H
@@ -0,0 +1,78 @@
+/*------------------------------- 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 __geometric_H__
+#define __geometric_H__ 
+
+#include "vtkFile.H"
+#include "triSurface.H"
+#include "multiTriSurface.H"
+#include "IOobject.H"
+
+namespace pFlow
+{
+
+
+template<typename ObjType>
+bool geomObjectToVTK(IOfileHeader& header, real time, fileSystem destPath, word bName)
+{
+
+	if( ObjType::TYPENAME() != header.objectType() )return false;
+
+	auto ioObjPtr = IOobject::make<ObjType>(header);
+
+	auto& data = ioObjPtr().template getObject<ObjType>();
+
+	vtkFile vtk(destPath, bName, time);
+
+	if(!vtk) return false;
+
+	Report(1)<<"Converting geometry to vtk."<<endReport;
+
+	if(!dataToVTK(vtk, data) )
+	{
+		fatalErrorInFunction<<
+		"  error in writing object "<<ioObjPtr().typeName() << " to folder " << destPath<<endl;
+		fatalExit;
+	}
+	
+	return true;
+}
+
+template<typename Type>
+bool dataToVTK(vtkFile& vtk, const Type& dataEntity)
+{
+	fatalErrorInFunction<<
+	"not implemented function!";
+	fatalExit;
+	return false;
+}
+
+template<>
+bool dataToVTK( vtkFile& vtk, const triSurface& surface );
+
+template<>
+bool dataToVTK( vtkFile& vtk, const multiTriSurface& surface );
+
+
+}
+
+#endif //__geometric_H__
diff --git a/utilities/pFlowToVTK/pFlowToVTK.C b/utilities/pFlowToVTK/pFlowToVTK.C
new file mode 100755
index 00000000..7087589a
--- /dev/null
+++ b/utilities/pFlowToVTK/pFlowToVTK.C
@@ -0,0 +1,142 @@
+/*------------------------------- 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 "systemControl.H"
+#include "pointFieldToVTK.H"
+#include "triSurfaceFieldToVTK.H"
+#include "timeFolder.H"
+#include "commandLine.H"
+#include "ranges.H"
+
+using pFlow::word;
+using pFlow::wordVector;
+using pFlow::geometryFolder__;
+using pFlow::timeFolder;
+using pFlow::fileSystem;
+using pFlow::wordList;
+using pFlow::IOfileHeader;
+using pFlow::objectFile;
+using pFlow::output;
+using pFlow::endl;
+using pFlow::multiTriSurface;
+using pFlow::commandLine;
+using pFlow::realCombinedRange;
+
+int main(int argc, char** argv )
+{
+	word outFolder = (pFlow::CWD()/word("VTK")).wordPath();
+
+	commandLine cmds(
+		"pFlowToVTK",
+		"Convrtes the saved pointField and geometry"
+		" date in time folders into vtk file format.");
+
+	wordVector times;
+		
+	bool noGoem = false;
+	cmds.add_flag(
+		"--no-geometry",
+		noGoem,
+		"Do not convert geometry to VTK file");
+	
+	bool noParticle = false;
+	cmds.add_flag("--no-particles",
+		noParticle,
+		"Do not convert particle fields to VTK file");
+
+	cmds.addOption("-o,--out-folder",
+		outFolder,
+		"path to output folder of VTK",
+		"path");
+
+	cmds.addOption(
+		"-t,--time",
+		times.vectorField(),
+		"a space separated lits of time folders, or a strided range begin:stride:end, or an interval begin:end",
+		" ");
+
+	if(!cmds.parse(argc, argv)) return 0;
+	
+	
+// this should be palced in each main 
+#include "initialize_Control.H"
+
+
+	timeFolder folders(Control);
+	fileSystem destFolder = fileSystem(outFolder)/geometryFolder__;
+	fileSystem destFolderField = fileSystem(outFolder);
+	wordList geomfiles{"triSurface"};
+
+	realCombinedRange validRange;
+	if( cmds.count("--time") )
+	{
+		if(!validRange.addRanges(times))
+		{
+			fatalExit;
+		}
+	}
+	else
+	{
+		validRange.addIntervalRange(folders.startTime(), folders.endTime());
+	}
+
+	do
+	{
+
+		if( !validRange.isMember( folders.time() ) )continue;
+		
+		output<< "time: " << cyanText( folders.time() )<<" s" <<endl;
+		if(!noGoem)
+		{	
+			fileSystem geomFolder = folders.folder()/geometryFolder__;
+			if(!pFlow::TSFtoVTK::convertTimeFolderTriSurfaceFields(geomFolder, folders.time(), destFolder, "surface"))
+			{
+				fatalExit;
+				return 1;
+			}	
+		}
+
+		if(!noParticle)
+		{
+			if( !pFlow::PFtoVTK::convertTimeFolderPointFields(
+				folders.folder(),
+				folders.time(),
+				destFolderField,
+				"sphereFields" )
+			)
+			{
+				fatalExit;
+			}	
+		}
+		
+		output<<endl;
+
+	}
+	while( folders++ );
+
+
+	output<< "\nFinished successfully.\n";
+
+	
+// this should be palced in each main 
+#include "finalize.H"
+	return 0;
+}
diff --git a/utilities/pFlowToVTK/pointFieldToVTK.H b/utilities/pFlowToVTK/pointFieldToVTK.H
new file mode 100755
index 00000000..826d2898
--- /dev/null
+++ b/utilities/pFlowToVTK/pointFieldToVTK.H
@@ -0,0 +1,361 @@
+/*------------------------------- 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 __pointFieldToVTK_H__
+#define __pointFieldToVTK_H__
+
+#include <regex>
+
+#include "vtkFile.H"
+#include "pointFields.H"
+#include "IOobject.H"
+
+namespace pFlow::PFtoVTK
+{
+
+template<typename IncludeMaskType>
+bool addInt64PointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	int64* 		field,
+	IncludeMaskType includeMask );
+
+template<typename IncludeMaskType>
+bool addRealPointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	real* 		field,
+	IncludeMaskType includeMask );
+
+template<typename IncludeMaskType>
+bool addRealx3PointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	realx3* 	field,
+	IncludeMaskType includeMask );
+
+bool regexCheck(word TYPENAME, word fieldType)
+{
+	std::regex match("pointField\\<([A-Za-z1-9_]*)\\,([A-Za-z1-9_]*)\\>");
+	std::smatch search1, search2;
+	if(!std::regex_match(fieldType, search1, match))return false;
+	if(!std::regex_match(TYPENAME, search2, match))return false;
+	if(search1.size()!=3)return false;
+	if(search1.size()!=search2.size())return false;
+	return search1[1] == search2[1];
+}
+
+template<typename Type>
+bool checkFieldType(word objectType)
+{
+	//if( pointField<VectorSingle,Type>::TYPENAME() == objectType )return true;
+	//if( pointField<VectorSingle,Type, HostSpace>::TYPENAME() == objectType ) return true;
+	//if( pointField<VectorDual, Type>::TYPENAME() == objectType )return true;
+	return regexCheck(pointField<VectorSingle,Type>::TYPENAME(), objectType);
+	
+}
+
+bool convertIntTypesPointField(
+	iOstream& os,
+	const IOfileHeader& header,
+	const pointStructure& pStruct )
+{
+	word objectType = header.objectType();
+
+	if( !(checkFieldType<int8>(objectType)  ||
+		checkFieldType<int16>(objectType) 	||
+		checkFieldType<int32>(objectType) 	||
+		checkFieldType<int64>(objectType) 	||
+		checkFieldType<uint32>(objectType) 	||
+		checkFieldType<label>(objectType))
+		)
+	{
+		return false;
+	}
+
+	auto objField = IOobject::make<int64PointField_H>
+	(
+		header,
+		pStruct,
+		static_cast<int64>(0)
+	);
+
+	auto& Field = objField().getObject<int64PointField_H>();
+
+	int64* data = Field.hostVectorAll().data();
+	
+	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk.\n";
+
+	return addInt64PointField(
+		os,
+		header.objectName(),
+		pStruct.numActive(),
+		data,
+		pStruct.activePointsMaskH() );
+
+}
+
+bool convertRealTypePointField(
+	iOstream& os,
+	const IOfileHeader& header,
+	const pointStructure& pStruct)
+{
+	word objectType = header.objectType();
+
+	if(!checkFieldType<real>(objectType))return false;
+
+	auto objField = IOobject::make<realPointField_H>
+	(
+		header,
+		pStruct,
+		static_cast<real>(0)
+	);
+
+	auto& Field = objField().getObject<realPointField_H>();
+
+	real* data = Field.hostVectorAll().data();
+	
+	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+
+	return addRealPointField(
+		os,
+		header.objectName(),
+		pStruct.numActive(),
+		data,
+		pStruct.activePointsMaskH() );
+}
+
+bool convertRealx3TypePointField(
+	iOstream& os,
+	const IOfileHeader& header,
+	const pointStructure& pStruct)
+{
+	word objectType = header.objectType();
+
+	if(!checkFieldType<realx3>(objectType))return false;
+
+	auto objField = IOobject::make<realx3PointField_H>
+	(
+		header,
+		pStruct,
+		static_cast<real>(0)
+	);
+
+	auto& Field = objField().getObject<realx3PointField_H>();
+
+	realx3* data = Field.hostVectorAll().data();
+	
+	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+
+	return addRealx3PointField(
+		os,
+		header.objectName(),
+		pStruct.numActive(),
+		data,
+		pStruct.activePointsMaskH() );
+}
+
+
+template<typename IncludeMaskType>
+bool addUndstrcuturedGridField(
+	iOstream& 	os,
+	int32 		numActivePoints,
+	realx3* 	position,
+	IncludeMaskType includeMask)
+{
+	
+	
+
+	auto [iFirst, iLast] = includeMask.activeRange();
+
+	os<< "DATASET UNSTRUCTURED_GRID\n";
+	os<< "POINTS "<< numActivePoints << " float\n"; 
+	
+	if(numActivePoints==0) return true;
+	
+	for(int32 i=iFirst; i<iLast; ++i)
+	{
+		if(includeMask(i))
+			os<< position[i].x()<<' '<< position[i].y()<<' '<<position[i].z()<<'\n';
+	}
+
+	os<<"CELLS "<< numActivePoints<<' '<< 2*numActivePoints<<'\n';
+	for(int32 i=0; i<numActivePoints; i++)
+	{
+		os<< 1 <<' '<< i<<'\n';
+	}
+
+	os<<"CELL_TYPES "<< numActivePoints<<'\n';
+
+	for(int32 i=0; i<numActivePoints; i++)
+	{
+		os<< 1 <<'\n';
+	}
+	
+	os << "POINT_DATA " << numActivePoints << endl;
+
+	return true;
+}
+
+
+template<typename IncludeMaskType>
+bool addInt64PointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	int64* 		field,
+	IncludeMaskType includeMask )
+{
+	if(numActivePoints==0) return true;
+
+	auto [iFirst, iLast] = includeMask.activeRange();
+
+	os << "FIELD FieldData 1\n"<<
+	fieldName << " 1 " << numActivePoints << " int\n";
+	for(int32 i=iFirst; i<iLast; ++i)
+	{
+		if(includeMask(i))
+			os<< field[i] <<'\n';
+	}
+
+	return true;
+}
+
+template<typename IncludeMaskType>
+bool addRealPointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	real* 		field,
+	IncludeMaskType includeMask )
+{
+	if(numActivePoints==0) return true;
+
+	auto [iFirst, iLast] = includeMask.activeRange();
+
+	os << "FIELD FieldData 1\n"<<
+	fieldName << " 1 " << numActivePoints << " float\n";
+	for(int32 i=iFirst; i<iLast; ++i)
+	{
+		if(includeMask(i))
+			os<< field[i] <<'\n';
+	}
+	return true;
+}
+
+template<typename IncludeMaskType>
+bool addRealx3PointField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		numActivePoints,
+	realx3* 	field,
+	IncludeMaskType includeMask )
+{
+	if(numActivePoints==0) return true;
+
+	auto [iFirst, iLast] = includeMask.activeRange();
+
+	os << "FIELD FieldData 1\n"<<
+	fieldName << " 3 " << numActivePoints << " float\n";
+	for(int32 i=iFirst; i<iLast; ++i)
+	{
+		if(includeMask(i))
+			os<< field[i].x()<<' '<< field[i].y()<<' '<<field[i].z()<<'\n';
+	}
+
+	return true;
+}
+
+
+bool convertTimeFolderPointFields( 
+	fileSystem timeFolder,
+	real time,
+	fileSystem destPath,
+	word bName)
+{
+
+	// check if pointStructure exist in this folder 
+	IOfileHeader pStructHeader(
+		objectFile(
+			pointStructureFile__,
+			timeFolder,
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS)
+		);
+
+	if( !pStructHeader.headerOk(true) )
+	{
+		output<<yellowColor<<"Time folder "<< timeFolder <<
+		" does not contain any pStructure data file. Skipping this folder . . ."
+		<<defaultColor<<nl;
+		return true;
+	}
+
+	vtkFile vtk(destPath, bName, time);
+
+	if(!vtk) return false;
+
+	auto pStructObjPtr = IOobject::make<pointStructure>(pStructHeader);
+	auto& pStruct = pStructObjPtr().getObject<pointStructure>();
+
+	// get a list of files in this timeFolder;
+
+	auto posVec = std::as_const(pStruct).pointPosition().hostVectorAll();
+	auto* pos = posVec.data();
+
+	Report(1)<<"Writing pointStructure to vtk file."<<endReport;
+	addUndstrcuturedGridField(
+		vtk(),
+		pStruct.numActive(),
+		pos,
+		pStruct.activePointsMaskH());
+
+	auto fileList = containingFiles(timeFolder);
+
+
+	for(auto& file:fileList)
+	{
+		IOfileHeader fieldHeader(
+			objectFile(
+			file.wordPath(),
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS) );
+
+		if( fieldHeader.headerOk(true) )
+		{
+			convertIntTypesPointField(vtk(), fieldHeader, pStruct);
+			convertRealTypePointField(vtk(), fieldHeader, pStruct);
+			convertRealx3TypePointField(vtk(), fieldHeader, pStruct);
+		}
+	}
+
+	return true;
+}
+ 
+}
+
+
+
+#endif 
diff --git a/utilities/pFlowToVTK/triSurfaceFieldToVTK.H b/utilities/pFlowToVTK/triSurfaceFieldToVTK.H
new file mode 100644
index 00000000..689435d1
--- /dev/null
+++ b/utilities/pFlowToVTK/triSurfaceFieldToVTK.H
@@ -0,0 +1,245 @@
+/*------------------------------- 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 __triSurfaceFieldToVTK_H__
+#define __triSurfaceFieldToVTK_H__
+
+#include <regex>
+
+#include "vtkFile.H"
+#include "triSurface.H"
+#include "multiTriSurface.H"
+#include "triSurfaceFields.H"
+#include "IOobject.H"
+
+namespace pFlow::TSFtoVTK
+{
+
+bool regexCheck(word TYPENAME, word fieldType)
+{
+	std::regex match("triSurfaceField\\<([A-Za-z1-9_]*)\\,([A-Za-z1-9_]*)\\>");
+	std::smatch search1, search2;
+	if(!std::regex_match(fieldType, search1, match))return false;
+	if(!std::regex_match(TYPENAME, search2, match))return false;
+	if(search1.size()!=3)return false;
+	if(search1.size()!=search2.size())return false;
+	return search1[1] == search2[1];
+}
+
+template<typename Type>
+bool checkFieldType(word objectType)
+{
+	//if( pointField<VectorSingle,Type>::TYPENAME() == objectType )return true;
+	//if( pointField<VectorSingle,Type, HostSpace>::TYPENAME() == objectType ) return true;
+	//if( pointField<VectorDual, Type>::TYPENAME() == objectType )return true;
+	return regexCheck(triSurfaceField<VectorSingle,Type>::TYPENAME(), objectType);
+	
+}
+
+template<typename Type>
+bool triDataToVTK(iOstream& os, const Type& dataEntity)
+{
+	fatalErrorInFunction<<
+	"not implemented function!";
+	fatalExit;
+	return false;
+}
+
+template<>
+bool triDataToVTK(iOstream& os, const triSurface& surface )
+{
+	auto nP = surface.numPoints();
+	auto hPoints = surface.points().hostVector();
+
+	os << "DATASET POLYDATA" << endl;
+	os << "POINTS " << nP << " float" << endl;
+
+
+	for ( auto i=0; i<nP; i++ )
+	{
+		os << hPoints[i].x()  << " " << hPoints[i].y() << " " << hPoints[i].z() << endl;
+	}
+
+	auto nV = surface.numTriangles();
+	auto hVertices = surface.vertices().hostVector();
+
+	os << "POLYGONS " << nV << "  " << 4*nV << endl;
+
+	for(auto i=0; i<nV; i++)
+	{
+		os<< 3 <<" "<< hVertices[i].x() << " " << hVertices[i].y() <<" "<<hVertices[i].z()<<endl;
+	}
+
+	return true;
+}
+
+template<>
+bool triDataToVTK(iOstream& os, const multiTriSurface& surface )
+{
+	auto nP = surface.numPoints();
+	auto hPoints = surface.points().hostVector();
+
+	os << "DATASET UNSTRUCTURED_GRID" << endl;
+	os << "POINTS " << nP << " float" << endl;
+
+
+	for ( auto i=0; i<nP; i++ )
+	{
+		os << hPoints[i].x()  << " " << hPoints[i].y() << " " << hPoints[i].z() << endl;
+	}
+
+	auto nV = surface.numTriangles();
+	auto hVertices = surface.vertices().hostVector();
+
+	os<<"CELLS "<< nV<<' '<< 4*nV<<'\n';
+	//os << "POLYGONS " << nV << "  " << 4*nV << endl;
+
+	for(auto i=0; i<nV; i++)
+	{
+		os<< 3 <<" "<< hVertices[i].x() << " " << hVertices[i].y() <<" "<<hVertices[i].z()<<endl;
+	}
+
+	os<<"CELL_TYPES "<< nV<<'\n';
+
+	for(int32 i=0; i<nV; i++)
+	{
+		os<< 5 <<'\n';
+	}
+	
+	os << "CELL_DATA " << nV << endl;
+
+	return true;
+}
+
+bool addRealx3TriSurfaceField(
+	iOstream& 	os,
+	word 		fieldName,
+	int32 		size,
+	realx3* 	field )
+{
+	if(size==0) return true;
+
+
+	os << "FIELD FieldData 1\n"<<
+	fieldName << " 3 " << size << " float\n";
+	for(int32 i=0; i<size; ++i)
+	{
+			os<< field[i].x()<<' '<< field[i].y()<<' '<<field[i].z()<<'\n';
+	}
+
+	return true;
+}
+
+bool convertRealx3TypetriSurfaceField(
+	iOstream& os,
+	const IOfileHeader& header,
+	const multiTriSurface& tSurface)
+{
+	word objectType = header.objectType();
+
+	if(!checkFieldType<realx3>(objectType))return false;
+
+	auto objField = IOobject::make<realx3TriSurfaceField_H>
+	(
+		header,
+		tSurface,
+		static_cast<real>(0)
+	);
+
+	auto& Field = objField().getObject<realx3TriSurfaceField_H>();
+
+	realx3* data = Field.hostVectorAll().data();
+	
+	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+
+	return addRealx3TriSurfaceField(
+		os,
+		header.objectName(),
+		tSurface.size(),
+		data );
+}
+
+bool convertTimeFolderTriSurfaceFields( 
+	fileSystem timeFolder,
+	real time,
+	fileSystem destPath,
+	word bName)
+{
+
+	// check if pointStructure exist in this folder 
+	IOfileHeader triSurfaeHeader(
+		objectFile(
+			triSurfaceFile__,
+			timeFolder,
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS)
+		);
+
+	if( !triSurfaeHeader.headerOk(true) )
+	{
+		output<<yellowText("Time folder "<< timeFolder <<
+		" does not contain any triSurface data file. Skipping this folder . . ."
+		)<<nl;
+		return true;
+	}
+
+	vtkFile vtk(destPath, bName, time);
+
+	if(!vtk) return false;
+
+	auto triSurfaceObjPtr = IOobject::make<multiTriSurface>(triSurfaeHeader);
+	auto& tSurface = triSurfaceObjPtr().getObject<multiTriSurface>();
+
+	// get a list of files in this timeFolder;
+	Report(1)<<"Wrting triSurface mesh/Geometry to vtk file."<<endReport;
+	if(!triDataToVTK(vtk(), tSurface))
+	{
+		fatalErrorInFunction<<
+		"error in writing triSurface data to vtk file "<< vtk.fileName()<<endl;
+		return false;
+	}
+
+
+	auto fileList = containingFiles(timeFolder);
+
+
+	for(auto& file:fileList)
+	{
+		IOfileHeader fieldHeader(
+			objectFile(
+			file.wordPath(),
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS) );
+
+		if( fieldHeader.headerOk(true) )
+		{
+			//output<<"object file type is "<<fieldHeader.objectType()<<endl;
+			convertRealx3TypetriSurfaceField(vtk(), fieldHeader, tSurface);
+		}
+	}
+
+	return true;
+}
+
+}
+
+#endif 
diff --git a/utilities/pFlowToVTK/vtkFile.C b/utilities/pFlowToVTK/vtkFile.C
new file mode 100755
index 00000000..d6388a6b
--- /dev/null
+++ b/utilities/pFlowToVTK/vtkFile.C
@@ -0,0 +1,75 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#include "vtkFile.H"
+
+bool pFlow::vtkFile::openStream()
+{
+	oStream_ = makeUnique<oFstream>( fileName() );
+	if( !oStream_ )return false;
+	return writeHeader();
+}
+
+bool pFlow::vtkFile::vtkFile::writeHeader()
+{
+
+	if(!oStream_) return false;
+
+	oStream_() << "# vtk DataFile Version 2.0" << endl;
+	
+	oStream_() << "vtk file for time : " << time_ << endl;
+
+	oStream_() << "ASCII" << endl;
+
+	if( oStream_().fail() ) return false;
+
+	return true;
+
+}
+
+pFlow::vtkFile::vtkFile
+(
+	const fileSystem dir,
+	const word& bName,
+	real time
+)
+:
+	dirPath_(dir),
+	baseName_(bName),
+	time_(time)
+{
+
+	if(!openStream())
+	{
+		fatalErrorInFunction <<
+		"  error in creating vtkFile "<<fileName()<<endl;
+		fatalExit;
+	}
+
+}
+
+
+pFlow::fileSystem pFlow::vtkFile::fileName()const
+{
+	word fName = baseName_ +"-" + int322Word(10000*time_) + ".vtk";
+	return dirPath_ +fName;
+}
+
diff --git a/utilities/pFlowToVTK/vtkFile.H b/utilities/pFlowToVTK/vtkFile.H
new file mode 100755
index 00000000..0d5861bf
--- /dev/null
+++ b/utilities/pFlowToVTK/vtkFile.H
@@ -0,0 +1,90 @@
+/*------------------------------- 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 __vtkFile_H__
+#define __vtkFile_H__ 
+
+#include "types.H"
+#include "uniquePtr.H"
+#include "fileSystem.H"
+#include "streams.H"
+
+namespace pFlow
+{
+
+
+class vtkFile
+{
+protected:
+
+	fileSystem dirPath_;
+
+	word       baseName_;
+
+	real  	   time_ = 0.0;
+
+	uniquePtr<oFstream>  oStream_= nullptr;
+
+	bool openStream();
+
+	virtual bool writeHeader();
+
+public:
+
+	vtkFile(const fileSystem dir, const word& bName, real time);
+
+	virtual ~vtkFile() = default;
+
+	inline oFstream& operator ()()
+	{
+		if(!oStream_)
+		{
+			if(!openStream())
+			{
+				fatalErrorInFunction<<
+				"  error in opening vtkFile "<< fileName() <<endl;
+				fatalExit;
+			}
+		}
+		return oStream_();
+	}
+
+	inline explicit operator bool() const
+	{
+		if(!oStream_)return false;
+		if(oStream_().fail())return false;
+		return true;
+	}
+
+	inline bool operator!()const
+	{
+		if( !oStream_) return true;
+		if( oStream_().fail() )return true;
+		
+		return false;
+	}
+
+	virtual fileSystem fileName()const;	
+
+};
+
+}
+
+#endif //__vtkFile_H__
diff --git a/utilities/particlesPhasicFlow/CMakeLists.txt b/utilities/particlesPhasicFlow/CMakeLists.txt
new file mode 100644
index 00000000..ee336506
--- /dev/null
+++ b/utilities/particlesPhasicFlow/CMakeLists.txt
@@ -0,0 +1,11 @@
+
+set(source_files
+particlesPhasicFlow.C
+positionParticles/positionParticles.C 
+positionOrdered/positionOrdered.C 
+positionRandom/positionRandom.C
+empty/empty.C
+)
+set(link_lib phasicFlow Kokkos::kokkos Interaction)
+
+pFlow_make_executable_install(particlesPhasicFlow source_files link_lib)
diff --git a/utilities/particlesPhasicFlow/empty/empty.C b/utilities/particlesPhasicFlow/empty/empty.C
new file mode 100755
index 00000000..ca2f8b1a
--- /dev/null
+++ b/utilities/particlesPhasicFlow/empty/empty.C
@@ -0,0 +1,36 @@
+/*------------------------------- 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 "empty.H"
+
+
+
+pFlow::empty::empty(
+	const dictionary& dict
+)
+:
+	positionParticles(dict),
+	position_
+	(
+		maxNumberOfParticles_, RESERVE()
+	)
+{
+	position_.clear();
+}
\ No newline at end of file
diff --git a/utilities/particlesPhasicFlow/empty/empty.H b/utilities/particlesPhasicFlow/empty/empty.H
new file mode 100755
index 00000000..3aaeef20
--- /dev/null
+++ b/utilities/particlesPhasicFlow/empty/empty.H
@@ -0,0 +1,94 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __empty_H__
+#define __empty_H__
+
+#include "positionParticles.H"
+#include "box.H"
+
+namespace pFlow
+{
+
+
+class empty 
+:
+	public positionParticles
+{
+protected:
+
+	dictionary 	emptyDict_;
+
+
+	realx3Vector 	position_;
+
+public:
+
+	// - type Info
+	TypeName("empty");
+
+	empty(const dictionary& dict);
+
+	// - add this class to vCtor selection table 
+	add_vCtor(
+		positionParticles,
+		empty,
+		dictionary);
+
+	virtual ~empty() = default;
+
+	//// - Methods 
+
+	virtual label numPoints()const
+	{
+		return 0;
+	}
+
+	virtual label size()const
+	{
+		return 0;
+	}
+
+	real maxDiameter() const override
+	{
+		return 1.0;
+	}
+
+	// - const access to position
+	virtual const realx3Vector& position()const 
+	{
+		return position_;
+	}
+
+	// - access to position 
+	virtual realx3Vector& position()
+	{
+		return position_;
+	}
+
+	
+};	
+
+
+}
+
+
+
+#endif // __empety_H__
diff --git a/utilities/particlesPhasicFlow/particlesPhasicFlow.C b/utilities/particlesPhasicFlow/particlesPhasicFlow.C
new file mode 100755
index 00000000..7d0a1f48
--- /dev/null
+++ b/utilities/particlesPhasicFlow/particlesPhasicFlow.C
@@ -0,0 +1,195 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#include "positionParticles.H"
+#include "pointStructure.H"
+#include "setFields.H"
+#include "systemControl.H"
+#include "commandLine.H"
+
+using pFlow::output;
+using pFlow::endl;
+using pFlow::dictionary;
+using pFlow::uniquePtr;
+using pFlow::IOobject;
+using pFlow::objectFile;
+using pFlow::fileSystem;
+using pFlow::pointStructure;
+using pFlow::pointStructureFile__;
+using pFlow::positionParticles;
+using pFlow::commandLine;
+using pFlow::setFieldList;
+
+int main( int argc, char* argv[] )
+{
+
+	commandLine cmds(
+		"createParticles",
+		"Read the dictionary createParticles and create particles"
+		" based on the two sub-dictionaries positionParticles and setFields.\n"
+		"First executespositionParticles and then setFields, except "
+		"otherwise selected in the command line.");
+
+			
+	bool positionOnly = false;
+	cmds.add_flag(
+		"--positionParticles-only",
+		positionOnly,
+		"Exectue the positionParticles part only and store the created "
+		"pointStructure in the time folder.");
+	
+	bool setOnly = false;
+	cmds.add_flag("--setFields-only",
+		setOnly,
+		"Exectue the setFields part only. Read the pointStructure from "
+		"time folder and setFields and save the result in the same time folder.");
+
+
+	if(!cmds.parse(argc, argv)) return 0;
+	
+	if(setOnly && positionOnly)
+	{
+		Err<<
+		"Options --positionParticles-only and --setFields-only cannot be used simeltanuously. \n"<<endErr;
+		return 1;
+	}	
+
+// this should be palced in each main 
+#include "initialize_Control.H"
+
+	auto objCPDict = IOobject::make<dictionary>
+	(
+		objectFile
+		(
+			"createParticles",
+			Control.settings().path(),
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
+		),
+		"createParticles",
+		true
+	);
+
+	auto& cpDict = objCPDict().getObject<dictionary>();
+
+	uniquePtr<IOobject> pStructObj{nullptr};
+
+	if(!setOnly)
+	{
+
+		// position particles based on the dict content 
+		Report(0)<< "Positioning points . . . \n"<<endReport;
+		auto pointPosition = positionParticles::create(cpDict.subDict("positionParticles"));
+
+		fileSystem pStructPath = Control.time().path()+pointStructureFile__;
+
+		auto finalPos = pointPosition().getFinalPosition();
+
+		
+		pStructObj = IOobject::make<pointStructure>
+		(
+			objectFile
+			(
+				pointStructureFile__,
+				Control.time().path(),
+				objectFile::READ_NEVER,
+				objectFile::WRITE_ALWAYS			
+			),
+			finalPos
+		);
+
+		auto& pSruct = pStructObj().getObject<pointStructure>();
+
+		Report(1)<< "Created pStruct with "<< pSruct.size() << " points and capacity "<<
+		pSruct.capacity()<<" . . ."<< endReport;
+
+		Report(1)<< "Writing pStruct to " << pStructObj().path() << endReport<<endl<<endl;
+
+		if( !pStructObj().write())
+		{
+			fatalErrorInFunction<<
+			"Error in writing to file. \n ";
+			return 1; 
+		}
+	}else
+	{
+		pStructObj = IOobject::make<pointStructure>
+		(
+			objectFile
+			(
+				pointStructureFile__,
+				Control.time().path(),
+				objectFile::READ_ALWAYS,
+				objectFile::WRITE_NEVER			
+			)
+		);
+	}
+
+	
+
+	if(!positionOnly)
+	{
+
+		auto& pStruct = pStructObj().getObject<pointStructure>();
+
+		auto& sfDict = cpDict.subDict("setFields");
+
+		setFieldList defValueList(sfDict.subDict("defaultValue"));
+
+		for(auto& sfEntry: defValueList)
+		{
+			if( !sfEntry.setPointFieldDefaultValueNewAll(Control.time(), pStruct, true))
+			{
+				Err<< "\n error occured in setting default value fields.\n"<<endErr;
+				return 1;
+			}
+		}
+		
+		output<<endl;
+
+		auto& selectorsDict = sfDict.subDict("selectors");
+
+		auto selNames = selectorsDict.dictionaryKeywords();
+
+		for(auto name: selNames)
+		{
+			Report(1)<< "Applying selector " << greenText(name) <<endReport;
+			
+			if(
+				!applySelector(Control, pStruct, selectorsDict.subDict(name))
+			 )
+			{
+				Err<<"\n error occured in setting selector. \n"<<endErr;
+				return 1;
+			}
+			output<<endl;
+		}
+	}
+
+	Control.time().write(true);
+	Report(0)<< greenText("\nFinished successfully.\n")<<endReport;
+
+
+// this should be palced in each main 
+#include "finalize.H"
+
+	return 0;
+} 
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.C b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.C
new file mode 100755
index 00000000..ea45cce8
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.C
@@ -0,0 +1,161 @@
+/*------------------------------- 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 "positionOrdered.H"
+#include "error.H"
+
+
+#include "streams.H"
+
+
+bool pFlow::positionOrdered::findAxisIndex()
+{
+	if( axisOrder_.size() != 3 )
+	{
+		fatalErrorInFunction <<
+		"axisOrder should have 3 components, but " << axisOrder_.size() << 
+		" has been provided. \n";
+		return false;
+	}
+
+
+	if( axisOrder_[0] == axisOrder_[1] || 
+		axisOrder_[0] == axisOrder_[2] ||
+		axisOrder_[1] == axisOrder_[2] )
+	{
+		fatalErrorInFunction<<
+		"invalid/repeated axis names in axisOrder. This is provided: " << axisOrder_ <<endl;
+		return false;
+	}
+
+	realx3 uV[3];
+	size_t i=0;
+	for(auto& ca: axisOrder_)
+	{
+		if(ca == "x")
+		{
+			uV[i] = realx3(1.0, 0.0, 0.0);
+		}
+		else if(ca == "y")
+		{
+			uV[i] = realx3(0.0, 1.0, 0.0);
+		}
+		else if(ca == "z")
+		{
+			uV[i] = realx3(0.0, 0.0, 1.0);
+		}
+		else
+		{
+			fatalErrorInFunction << 
+			"unknown name for axis in axisOrder: " << ca <<endl; 
+			return false;
+		}
+		i++;
+	}
+
+	uVector1_ = uV[0];
+	uVector2_ = uV[1];
+	uVector3_ = uV[2];
+
+	return true;
+}
+
+bool pFlow::positionOrdered::positionPointsOrdered()
+{
+	position_.clear();
+
+	realx3 dl(diameter_);
+	auto minP = static_cast<real>(0.5)*dl + box_.minPoint();
+	auto maxP = static_cast<real>(0.5)*dl + box_.maxPoint();
+
+	auto cntr = minP;
+
+	size_t n = 0;
+	while( n < numPoints_ )
+	{
+		position_.push_back(cntr);
+
+		cntr += dl*uVector1_;
+		
+		if( dot(uVector1_, cntr) > dot(uVector1_, maxP) )
+		{
+			cntr = (minP*uVector1_) + ( (cntr+dl) * uVector2_) + (cntr*uVector3_);
+
+			if( dot(uVector2_, cntr) > dot(uVector2_, maxP) )
+			{
+				cntr = (cntr*uVector1_) + (minP*uVector2_) + ((cntr+dl)*uVector3_);
+
+				if( dot(uVector3_,cntr) > dot(uVector3_, maxP) )
+				{
+					fatalErrorInFunction<<
+					"positioned " << n << " points in the domain and it is full. \n" <<
+					"request to position "<< numPoints_<< " points has failed.\n";
+					return false;
+				}
+			}
+		}
+		n++;
+	}
+
+	return true;
+}
+
+pFlow::positionOrdered::positionOrdered
+(
+	const dictionary& dict
+)
+:
+	positionParticles(dict),
+	poDict_
+	(  
+		dict.subDict("positionOrderedInfo")
+	),
+	diameter_
+	(
+		poDict_.getVal<real>("diameter")
+	),
+	numPoints_
+	(
+		poDict_.getVal<size_t>("numPoints")
+	),
+	axisOrder_
+	(
+		poDict_.getValOrSet("axisOrder", wordList{"x", "y", "z"})
+	),
+	box_
+	(
+		poDict_.subDict("box")
+	),
+	position_
+	(
+		maxNumberOfParticles_, RESERVE()
+	)
+{
+	
+	if( !findAxisIndex() )
+	{
+		fatalExit;
+	}
+
+	if(!positionPointsOrdered())
+	{
+		fatalExit;
+	}
+}
\ No newline at end of file
diff --git a/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.H b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.H
new file mode 100755
index 00000000..a2c10b81
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionOrdered/positionOrdered.H
@@ -0,0 +1,115 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __positionOrdered_H__
+#define __positionOrdered_H__
+
+#include "positionParticles.H"
+#include "box.H"
+
+namespace pFlow
+{
+
+
+class positionOrdered 
+:
+	public positionParticles
+{
+protected:
+
+	dictionary 	poDict_;
+
+	real 		diameter_;
+
+	size_t 		numPoints_;
+	
+	wordList 	axisOrder_;
+
+	box 		box_;
+
+	// - unit vector of the first axis
+	realx3 		uVector1_;
+
+	// - unit vector of the second axis
+	realx3 		uVector2_;
+
+	// - unit vector of the third axis
+	realx3 		uVector3_;
+
+
+	realx3Vector 	position_;
+
+	bool findAxisIndex();
+
+	bool positionPointsOrdered();
+
+public:
+
+	// - type Info
+	TypeName("positionOrdered");
+
+	positionOrdered(const dictionary& dict);
+
+	// - add this class to vCtor selection table 
+	add_vCtor(
+		positionParticles,
+		positionOrdered,
+		dictionary);
+
+	virtual ~positionOrdered() = default;
+
+	//// - Methods 
+
+	virtual label numPoints()const
+	{
+		return position_.size();
+	}
+
+	virtual label size()const
+	{
+		return position_.size();
+	}
+
+	real maxDiameter() const override
+	{
+		return diameter_;
+	}
+	
+	// - const access to position
+	virtual const realx3Vector& position()const 
+	{
+		return position_;
+	}
+
+	// - access to position 
+	virtual realx3Vector& position()
+	{
+		return position_;
+	}
+
+	
+};	
+
+
+}
+
+
+
+#endif // __positionOrdered_H__
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.C b/utilities/particlesPhasicFlow/positionParticles/positionParticles.C
new file mode 100755
index 00000000..d812b670
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.C
@@ -0,0 +1,110 @@
+/*------------------------------- 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 "positionParticles.H"
+#include "cells.H"
+#include "contactSearchFunctions.H"
+
+#include "streams.H"
+
+
+pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(realx3Vector& position)const
+{
+	struct indexMorton
+	{
+		size_t morton;
+		size_t index;
+	};
+
+	realx3 minP = min(position);
+	realx3 maxP = max(position);
+	real cellsize = maxDiameter();
+	cells<size_t> allCells( box(minP, maxP), cellsize);
+
+	Vector<indexMorton> indMor(position.size(),RESERVE());
+
+	indMor.clear();
+
+	size_t ind=0;
+	for(const auto& p:position)
+	{
+		auto cellInd = allCells.pointIndex(p);
+		indMor.push_back( 
+			{ xyzToMortonCode64(cellInd.x(), cellInd.y(), cellInd.z()),
+			ind++});
+	}
+
+	Info<<"Performing morton sorting."<<endInfo;
+	std::sort(
+		indMor.begin(),
+		indMor.end(),
+		[]( const indexMorton &lhs, const indexMorton &rhs){
+			return lhs.morton < rhs.morton; } );
+
+	realx3Vector sortedPos(position.capacity(), RESERVE());
+	sortedPos.clear();
+
+
+	for(auto& ind:indMor)
+	{
+		sortedPos.push_back( position[ind.index] );
+	}
+
+	return sortedPos;
+}
+
+
+pFlow::positionParticles::positionParticles
+(
+	const dictionary& dict
+)
+{
+	maxNumberOfParticles_ = dict.getValOrSet("maxNumberOfParticles", static_cast<size_t>(10000));
+	
+	mortonSorting_ = dict.getValOrSet("mortonSorting", Logical("Yes"));	
+}
+
+
+
+
+pFlow::uniquePtr<pFlow::positionParticles> pFlow::positionParticles::create(const dictionary & dict)
+{
+
+	word method = dict.getVal<word>("method");
+
+
+	if( dictionaryvCtorSelector_.search(method) )
+	{
+		return dictionaryvCtorSelector_[method] (dict);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< method << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			dictionaryvCtorSelector_
+		);
+		fatalExit;
+	}
+
+	return nullptr;
+}
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.H b/utilities/particlesPhasicFlow/positionParticles/positionParticles.H
new file mode 100755
index 00000000..41f8f7d6
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.H
@@ -0,0 +1,98 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __positionParticles_H__
+#define __positionParticles_H__
+
+#include "virtualConstructor.H"
+#include "Vectors.H"
+#include "dictionary.H"
+
+namespace pFlow
+{
+
+
+class positionParticles
+{
+protected:
+
+	size_t  		maxNumberOfParticles_ = 10000;
+
+	Logical 		mortonSorting_;
+
+	static const size_t numReports_ = 40;
+
+	realx3Vector sortByMortonCode(realx3Vector& position)const;
+	
+public:
+
+	// - type Info
+	TypeName("positionParticles");
+
+	positionParticles(const dictionary& dict);
+
+	create_vCtor
+	(
+		positionParticles,
+		dictionary,
+		(const dictionary& dict),
+		(dict)
+	);
+
+	virtual ~positionParticles() = default;
+
+	//// - Methods 
+
+	virtual label numPoints()const = 0;
+
+	virtual label size()const = 0;
+
+	virtual real maxDiameter() const = 0;
+
+	// - const access to position
+	virtual const realx3Vector& position()const = 0;
+	
+
+	// - access to position 
+	virtual realx3Vector& position() = 0;
+	
+	virtual realx3Vector getFinalPosition()
+	{
+		if(mortonSorting_)
+		{
+			return sortByMortonCode(position());
+		}
+		else
+		{
+			return position();
+		}
+	}
+
+
+	static
+	uniquePtr<positionParticles> create(const dictionary & dict);
+};	
+
+
+}
+
+
+
+#endif // __positionParticles_H__
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.C b/utilities/particlesPhasicFlow/positionRandom/positionRandom.C
new file mode 100755
index 00000000..628a1903
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.C
@@ -0,0 +1,260 @@
+/*------------------------------- 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 "positionRandom.H"
+#include "uniformRandomReal.H"
+#include "VectorSingles.H"
+#include "VectorDuals.H"
+#include "NBS.H"
+#include "unsortedPairs.H"
+
+
+
+namespace pFlow
+{
+
+using SearchType = NBS<DefaultExecutionSpace, int32> ;
+using ContainerType = unsortedPairs<DefaultExecutionSpace, int32>;
+
+
+void fillPoints(
+		uint numPoints,
+		realx3 minP,
+		realx3 maxP,
+		realx3Vector_HD& points,
+		int32Vector_HD& flags );
+
+int32 findCollisions(
+		ContainerType& pairs,
+		int32Vector_HD& flags);
+
+int32 findCollisions(int32 num, realx3Vector_HD& points, real diam)
+{
+	int32 res =0;
+	for(auto i=0; i<num;i++)
+	{
+		for(auto j=i+1; j<num; j++)
+		{
+			if(sphereSphereCheck(points[i],points[j],diam,diam))res++;
+		}
+	}
+
+	return res;
+}
+
+}
+
+
+bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
+{
+
+	realVector_D diameter(startNum , diameter_);
+	int32Vector_HD  flagHD(startNum, 0);
+	realx3Vector_HD positionHD(startNum);
+	
+	auto minP = box_.minPoint() + static_cast<real>(0.5)*realx3(diameter_);
+	auto maxP = box_.maxPoint() - static_cast<real>(0.5)*realx3(diameter_);
+
+	SearchType search(
+					box_,
+					diameter_,
+					positionHD.deviceVectorAll(),
+					diameter.deviceVectorAll());
+
+	ContainerType pairs(3*startNum);
+		
+	Report(1)<< "Positioning "<<
+		greenText("(Pass #"<< pass+1<<")")<< 
+		": started with "<< startNum <<" points."<<endReport;
+	
+	fillPoints(startNum, minP, maxP, positionHD, flagHD);
+		
+	search.broadSearch(pairs, range(0, startNum), true);
+
+	
+	int32 numCollisions = findCollisions(pairs, flagHD);
+	
+	
+	Report(2)<< "Positioned " << cyanText(startNum - numCollisions) <<
+	" without collision \n"<<endReport; 
+
+	if(startNum-numCollisions >= numPoints_ )
+	{
+		
+		Report(1)<<"Selected "<< cyanText(numPoints_)<< " for the final field.\n"<<endReport;
+
+		positionHD.syncViews();
+		position_.clear();
+		int32 n=0;
+		for(int32 i=0; i<startNum; i++)
+		{
+			if(flagHD[i] == 0 )
+			{
+				position_.push_back( positionHD[i]);
+				n++;
+				if(n==numPoints_)break;	
+			}
+			
+		}
+		
+		return true;
+	}
+	
+	
+	return false;
+	
+	
+}
+
+bool pFlow::positionRandom::positionPointsRandom()
+{
+	
+	position_.clear();
+
+	if(numPoints_ == 0)return true;
+	
+	size_t pass = 0;
+	int32 startNum = numPoints_;
+
+	while ( pass <maxIterations_)
+	{
+		if( positionOnePass(pass, startNum) )return true;
+		startNum = 1.1*startNum+1;
+		pass++;
+	}	
+	
+	
+	fatalErrorInFunction<< 
+	"  cannot position "<< numPoints_ << " in the domain in " << maxIterations_ << " iterations.\n" <<
+	"  you may increase maxIterations for positioning points.\n";
+
+	return false;
+}
+
+bool pFlow::positionRandom::inCollision
+(
+	const realx3 &cntr,
+	real diam
+)
+{
+	for(const auto& cp: position_)
+	{
+		if( length(cp-cntr) <= diam ) return true;
+	}
+
+	return false;
+}
+
+pFlow::positionRandom::positionRandom
+(
+	const dictionary& dict
+)
+:
+	positionParticles(dict),
+	prDict_
+	(  
+		dict.subDict("positionRandomInfo")
+	),
+	diameter_
+	(
+		prDict_.getVal<real>("diameter")
+	),
+	numPoints_
+	(
+		prDict_.getVal<size_t>("numPoints")
+	),
+	box_
+	(
+		prDict_.subDict("box")
+	),
+	maxIterations_
+	(
+		prDict_.getValOrSet("maxIterations", 10)
+	),
+	position_
+	(
+		maxNumberOfParticles_, RESERVE()
+	)
+{
+
+	reportInterval_ = max(numPoints_/numReports_, static_cast<size_t>(2));
+
+	if( !positionPointsRandom() )
+	{
+		fatalExit;
+	}
+}
+
+void pFlow::fillPoints(
+	uint numPoints,
+	realx3 minP,
+	realx3 maxP,
+	realx3Vector_HD& points,
+	int32Vector_HD& flags )
+{
+
+	uniformRandomReal rand;
+	
+
+	for(size_t i=0; i<numPoints; i++)
+	{
+		if(flags[i] == 0)
+		{
+			points[i] =rand(minP, maxP);	
+		}
+	}
+	points.modifyOnHost();
+	points.syncViews();
+}
+
+pFlow::int32 pFlow::findCollisions(
+		ContainerType& pairs,
+		int32Vector_HD& flags)
+{
+	auto allPairs = pairs.getPairs();
+	auto num = pairs.capacity();
+	auto dFlags = flags.deviceVector();
+	
+
+	int32 numCollisions = 0;
+
+	Kokkos::parallel_reduce(
+		"positionRandom::findCollisions",
+		num,
+		LAMBDA_HD(int32 i, int32& valueToUpdate){
+			if(allPairs.isValid(i))
+			{
+				auto pair = allPairs.getPair(i);
+				if( dFlags[pair.first] ==0 )
+				{
+					dFlags[pair.first] = 1;
+					valueToUpdate++;	
+				}
+			}
+		}, numCollisions);
+
+	flags.modifyOnDevice();
+	flags.syncViews();	
+
+	return numCollisions;
+}
+
+
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.H b/utilities/particlesPhasicFlow/positionRandom/positionRandom.H
new file mode 100755
index 00000000..1d145b5a
--- /dev/null
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.H
@@ -0,0 +1,113 @@
+/*------------------------------- 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 __positionOrdered_H__
+#define __positionOrdered_H__
+
+#include "positionParticles.H"
+#include "VectorSingles.H"
+
+#include "box.H"
+
+namespace pFlow
+{
+
+class positionRandom
+:
+	public positionParticles
+{
+protected:
+
+	dictionary 	prDict_;
+
+	real 		diameter_;
+
+	size_t 		numPoints_;
+	
+	box 		box_;
+
+	size_t 		maxIterations_;
+	
+	realx3Vector 	position_;
+
+	size_t 		reportInterval_;
+
+	bool positionOnePass(int32 pass, int32 startNum);
+	
+	bool positionPointsRandom();
+
+	bool inCollision(const realx3 &cntr, real diam);
+
+public:
+
+	// - type Info
+	TypeName("positionRandom");
+
+	positionRandom(const dictionary& dict);
+
+	// - add this class to vCtor selection table 
+	add_vCtor(
+		positionParticles,
+		positionRandom,
+		dictionary);
+
+	virtual ~positionRandom() = default;
+
+	//// - Methods 
+
+	virtual label numPoints()const
+	{
+		return position_.size();
+	}
+
+	virtual label size()const
+	{
+		return position_.size();
+	}
+
+	real maxDiameter() const override
+	{
+		return diameter_;
+	}
+
+	// - const access to position
+	virtual const realx3Vector& position()const 
+	{
+		return position_;
+	}
+
+	// - access to position 
+	virtual realx3Vector& position()
+	{
+		return position_;
+	}
+
+
+
+	
+};	
+
+
+}
+
+
+
+#endif // __positionOrdered_H__
diff --git a/utilities/particlesPhasicFlow/setFields.H b/utilities/particlesPhasicFlow/setFields.H
new file mode 100755
index 00000000..b56b4263
--- /dev/null
+++ b/utilities/particlesPhasicFlow/setFields.H
@@ -0,0 +1,56 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+
+#ifndef __setFields_H__
+#define __setFields_H__
+
+#include "pStructSelector.H"
+#include "pointFields.H"
+#include "setFieldList.H"
+#include "systemControl.H"
+
+namespace pFlow
+{
+
+
+bool applySelector(systemControl& control, const pointStructure& pStruct, const dictionary& selDict)
+{
+
+	
+	auto selector = pStructSelector::create(pStruct, selDict);
+
+	auto& selected = selector().selectedPoinsts();
+
+	int32IndexContainer selIndex(selected.data(), selected.size());
+
+	setFieldList sfList(selDict.subDict("fieldValue"));
+
+	for(auto& sfEntry:sfList)
+	{
+		if(!sfEntry.setPointFieldSelectedAll(control.time(), selIndex, true ))return false;
+	}
+	
+	return true;
+}
+
+} // pFlow
+
+#endif //__setFields_H__