diff --git a/src/phasicFlow/typeSelection/typeInfo.H b/src/phasicFlow/typeSelection/typeInfo.H
new file mode 100644
index 00000000..6f5596cf
--- /dev/null
+++ b/src/phasicFlow/typeSelection/typeInfo.H
@@ -0,0 +1,172 @@
+/*------------------------------- 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 __typeInfo_H__
+#define __typeInfo_H__
+
+#include <typeinfo>
+#include <cxxabi.h>
+
+#include "bTypes.H"
+
+
+#define has_static_member(name)		 								\
+	template<typename, typename> 									\
+	struct has_static_member_##name;								\
+																	\
+	template<typename testType, typename Ret, typename... Args>		\
+	struct has_static_member_##name<testType, Ret(Args...)> {		\
+    	template<typename U, U> struct Check; 						\
+																	\
+    template<typename U> 											\
+    static std::true_type Test(Check<Ret(*)(Args...), &U::name>*); 	\
+    template<typename U>											\
+    static std::false_type Test(...);								\
+    static const bool value = decltype(Test<testType>(0))::value;	\
+};
+
+
+
+#define TypeName(tName)  \
+	inline static word TYPENAME() {return tName; } \
+	virtual word typeName() const {return TYPENAME();} 
+
+#define TypeNameNV(tName)  \
+	inline static word TYPENAME() {return tName; } \
+	word typeName() const {return TYPENAME();} 
+
+
+#define TypeNameTemplate(tName, Type)										\
+	has_static_member(TYPENAME);											\
+	inline static word TYPENAME()											\
+	{ 																		\
+	if constexpr( has_static_member_TYPENAME<Type,word(void)>::value) 		\
+		{ return word(tName)+"<"+Type::TYPENAME()+">";} 						\
+	else																	\
+		return word(tName)+"<"+basicTypeName<Type>()+">";						\
+	}																		\
+	virtual word typeName() const { return TYPENAME();}
+
+#define TypeNameTemplate2(tName, Type1, Type2)								\
+	has_static_member(TYPENAME);											\
+	inline static word TYPENAME()											\
+	{ 																		\
+	if constexpr( has_static_member_TYPENAME<Type1,word(void)>::value) 		\
+		{ return word(tName)+"<"+Type1::TYPENAME()+","+Type2::TYPENAME()+">";}	\
+	else																	\
+		return word(tName)+"<"+basicTypeName<Type1>()+","+Type2::TYPENAME()+">";\
+	}																		\
+	virtual word typeName() const { return TYPENAME();}
+
+#define TypeNameTemplate3(tName, Type1, Type2, Type3)						\
+	inline static word TYPENAME()											\
+	{ 																		\
+		return word(tName)+"<"+Type1::TYPENAME()+","+Type2::TYPENAME()+","+Type3::TYPENAME()+">";\
+	}																		\
+	virtual word typeName() const { return TYPENAME();}
+
+// this is the non-virtual version
+#define TypeNameTemplateNV(tName, Type)										\
+	has_static_member(TYPENAME);											\
+	inline static word TYPENAME()											\
+	{ 																		\
+	if constexpr( has_static_member_TYPENAME<Type,word(void)>::value) 		\
+		{ return word(tName)+"<"+Type::TYPENAME()+">";} 					\
+	else																	\
+		return word(tName)+"<"+basicTypeName<Type>()+">";					\
+	return "noTYPE"; 														\
+	}																		\
+	inline word typeName() const { return TYPENAME();}
+
+
+#define TypeNameTemplateNV2(tName, Type, tName2)							\
+	has_static_member(TYPENAME);											\
+	inline static word TYPENAME()											\
+	{ 																		\
+	if constexpr( has_static_member_TYPENAME<Type,word(void)>::value) 		\
+		{ return word(tName)+"<"+Type::TYPENAME()+","+word(tName2)+">";} 	\
+	else																	\
+		return word(tName)+"<"+basicTypeName<Type>()+","+word(tName2)+">";	\
+	return "noTYPE"; 														\
+	}																		\
+	inline word typeName() const { return TYPENAME();}
+
+
+
+
+namespace pFlow
+{
+
+	template <typename T>
+	inline word basicTypeName()
+	{
+		int status;
+		auto& ti = typeid(T);
+		char* realname = abi::__cxa_demangle(ti.name(), 0, 0, &status);
+	  	word name(realname);
+	  	free(realname);
+	  	return name;
+	}
+
+	template<>
+	inline word basicTypeName<word>(){ return "word"; }
+
+	template<>
+	inline word basicTypeName<int64>(){ return "int64"; }
+
+	template<>
+	inline word basicTypeName<int32>(){ return "int32"; }
+
+	template<>
+	inline word basicTypeName<int16>(){ return "int16"; }
+
+	template<>
+	inline word basicTypeName<int8>(){ return "int8"; }
+
+	template<>
+	inline word basicTypeName<label>(){ return "label"; }
+
+	template<>
+	inline word basicTypeName<uint32>(){ return "uint32"; }
+
+	template<>
+	inline word basicTypeName<real>(){ return "real"; }
+
+		
+	// compare the overriden typeName of object with concrete TYPENAME 
+	// of Type1
+	template <typename Type1, typename Type2>
+	bool checkType(Type2* object)
+	{
+		return word(Type1::TYPENAME()) == object->typeName();
+	}
+
+	template <typename Type1, typename Type2>
+	bool checkType(Type2& object)
+	{
+		return word(Type1::TYPENAME()) == object.typeName();
+	}
+
+
+
+} // pFlow
+
+#endif
diff --git a/src/phasicFlow/typeSelection/virtualConstructor.H b/src/phasicFlow/typeSelection/virtualConstructor.H
new file mode 100644
index 00000000..5926c4a3
--- /dev/null
+++ b/src/phasicFlow/typeSelection/virtualConstructor.H
@@ -0,0 +1,82 @@
+/*------------------------------- 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 __virtualConstructor_H__
+#define __virtualConstructor_H__
+
+#include <functional>
+
+#include "types.H"
+#include "error.H"
+#include "Map.H"
+#include "uniquePtr.H"
+
+#define create_vCtor(baseClass,selectorName,argList,args) 								\
+																						\
+typedef std::function< uniquePtr<baseClass> argList > selectorName##FunctionType;		\
+typedef wordMap<selectorName##FunctionType>   selectorName##vCtorSelectorType;			\
+																						\
+inline static selectorName##vCtorSelectorType selectorName##vCtorSelector_; 			\
+																						\
+																						\
+template<typename dType>																\
+class create##selectorName##Callback 													\
+{																						\
+public: 																				\
+	create##selectorName##Callback ()		 											\
+	{																					\
+		auto success = 																	\
+			selectorName##vCtorSelector_.insertIf										\
+			(																			\
+					word(dType::TYPENAME()),											\
+					[&] argList -> uniquePtr<baseClass>									\
+					{ 																	\
+						return uniquePtr<baseClass>										\
+						( 																\
+							new dType args 												\
+						);																\
+					}																	\
+			);																			\
+																						\
+		if( !success )																	\
+		{																				\
+			fatalErrorInFunction														\
+			    << "Duplicate entry "<< dType::TYPENAME() 								\
+			    << " in virtual constructor table of "<< #baseClass						\
+			    << " with selector name " << #selectorName <<endl;						\
+		}																				\
+																						\
+	}																					\
+																						\
+	create##selectorName##Callback														\
+	(const create##selectorName##Callback&)= delete;									\
+	void operator=																		\
+	(const create##selectorName##Callback&)= delete;									\
+};
+
+
+
+#define add_vCtor(baseClass, derivedClass, selectorName) 	\
+															\
+inline static baseClass::create##selectorName##Callback<derivedClass> baseClass##derivedClass##selectorName##_;
+
+
+
+#endif // __virtualConstructor_H__
diff --git a/src/phasicFlow/types/basicTypes/Logical.C b/src/phasicFlow/types/basicTypes/Logical.C
new file mode 100755
index 00000000..29c450b1
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/Logical.C
@@ -0,0 +1,129 @@
+/*------------------------------- 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 "Logical.H"
+#include "error.H"
+#include "iIstream.H"
+#include "iOstream.H"
+
+
+pFlow::Logical::Logical(const word& l)
+{
+	if(!evaluteWord(l, s_, yesNoSet_ ))
+	{
+		fatalErrorInFunction<<
+		"  invalid input for Logical: "<< l << endl;
+		fatalExit;
+	}
+}
+
+pFlow::Logical::Logical(const char* ch)
+:
+	Logical(word(ch))
+{}
+
+bool pFlow::Logical::evaluteWord(const word& l, bool& b, int& yesNoSet )
+{
+	auto Ul = toUpper(l);
+	
+	for(int i=0; i<4; ++i)
+	{
+		if(toUpper(YesNo__[i][0]) == Ul)
+		{
+			b = true;
+			yesNoSet = i;
+			return true;
+		}
+		else if( toUpper(YesNo__[i][1]) == Ul )
+		{
+			b = false;
+			yesNoSet = i;
+			return true;
+		}
+	}
+	yesNoSet = 0;
+	return false;
+}
+
+bool pFlow::Logical::read(iIstream& is)
+{
+	token t(is);
+    word w;
+    if (!t.good())
+    {
+        ioErrorInFile(is.name(), is.lineNumber())
+            << "Bad token - could not get Logical value";
+        is.setBad();
+        return false;
+    }
+
+    if (t.isString())
+    {
+        w = t.stringToken(); 
+    }
+    else if(t.isWord() )
+    {
+        w = t.wordToken();
+    }
+    else
+    {
+        ioErrorInFile(is.name(), is.lineNumber())
+            << "Wrong token type - expected Logical value, found "
+            << t;
+        is.setBad();
+        return false;
+    }
+    
+    return evaluteWord(w, s_, yesNoSet_);
+}
+
+bool pFlow::Logical::write
+(
+	iOstream& os
+)const
+{
+	if(s_)
+	{
+		os<< YesNo__[yesNoSet_][0];
+	}
+	else
+	{
+		os<< YesNo__[yesNoSet_][1];
+	}
+	return os.check(FUNCTION_NAME);
+} 
+
+pFlow::iIstream& pFlow::operator>>( iIstream& is, Logical& L)
+{
+    if(!L.read(is))
+    {
+    	fatalExit;
+    }
+    return is;
+}
+
+pFlow::iOstream& pFlow::operator<< ( iOstream& os, const Logical& L)
+{
+	if(!L.write(os))
+	{
+		fatalExit;
+	}
+	return os;
+}
\ No newline at end of file
diff --git a/src/phasicFlow/types/basicTypes/Logical.H b/src/phasicFlow/types/basicTypes/Logical.H
new file mode 100755
index 00000000..22b8511b
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/Logical.H
@@ -0,0 +1,114 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __Logical_H__
+#define __Logical_H__
+
+#include "builtinTypes.H"
+#include "bTypesFunctions.H"
+#include "typeInfo.H"
+
+namespace pFlow
+{
+
+class iIstream;
+class iOstream;
+
+// allias for bool
+class Logical
+{
+protected:
+	bool s_ = false;
+
+	int  yesNoSet_ = 0;
+
+	inline static const word YesNo__[4][2] = {{"Yes", "No"},{"on","off"},{"true","false"}, {"Ok","No"}};
+
+	inline explicit Logical(bool s, int yns)
+	:	
+		s_(s),
+		yesNoSet_(yns)
+	{}
+
+public:
+
+	TypeNameNV("Logical");
+	
+
+	inline Logical(){}
+		
+
+	inline explicit Logical(bool s)
+	:
+	 	s_(s),
+	 	yesNoSet_(0)
+	{}
+
+	Logical(const word& l);
+	
+	Logical(const char* ch);
+
+	Logical(const Logical&) = default;
+
+	Logical(Logical&&) = default;
+
+	Logical& operator=(const Logical&) = default;
+
+	Logical& operator=(Logical&&) = default;
+
+	inline Logical& operator=(const bool& b)
+	{ 
+		s_ = b;
+		yesNoSet_ = 0;
+		return *this;
+	}
+	
+	inline bool operator()() const
+	{
+		return s_;
+	}
+	
+	inline explicit operator bool() const
+	{
+		return s_;
+	}
+
+	inline Logical operator!()const
+	{
+		return Logical(!s_, yesNoSet_);
+	}
+
+	//// IO operations
+	bool read(iIstream& is);
+
+	bool write(iOstream& os)const; 
+
+	bool static evaluteWord(const word& l, bool& b, int& yesNoSet );
+	
+
+};
+
+iIstream& operator>>( iIstream& is, Logical& L);
+
+iOstream& operator<<( iOstream& os, const Logical& L);
+
+} // pFlow
+
+#endif // __Logical_H__
diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.C b/src/phasicFlow/types/basicTypes/bTypesFunctions.C
new file mode 100755
index 00000000..a67fabf5
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.C
@@ -0,0 +1,320 @@
+/*------------------------------- 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 <algorithm>
+#include <sstream>
+#include <iomanip>
+
+#include "bTypesFunctions.H"
+
+
+pFlow::int32 pFlow::countChar(const word& s, const char c)
+{
+    return std::count(s.cbegin(), s.cend(), c);
+}
+
+
+pFlow::int32 pFlow::countChar( const char* s, const char c)
+{
+    return
+    (
+      s == nullptr
+      ? 0
+      : std::count(s, (s + std::char_traits<char>::length(s)), c)
+    );
+}
+
+pFlow::word pFlow::real2Fixed(const real & v, int32 numPrecision)
+{
+    std::stringstream ss;
+
+    ss << std::fixed << std::setprecision(numPrecision) << v;
+    return ss.str();
+}
+
+pFlow::word pFlow::real2Word(const real & v, int32 numPrecision)
+{
+    std::stringstream ss;
+    if( abs(v) < verySmallValue )
+    {
+        ss <<"0";
+    }
+    else
+    {
+        ss << std::setprecision(numPrecision) << v;
+    }
+    
+    return ss.str();
+}
+
+pFlow::word pFlow::int322Word(const int32 & v)
+{
+    std::stringstream ss;
+
+    ss << v;
+    return ss.str();
+}
+
+pFlow::word pFlow::toUpper(const word & inStr)
+{
+    word oStr(inStr);
+    transform(inStr.begin(), inStr.end(), oStr.begin(), ::toupper);
+    return oStr;
+}
+
+bool pFlow::isYes(const word & str)
+{
+    word s = toUpper(str);
+
+    if( s == "YES" || s == "OK" || s == "TRUE" || s == "ON" || s=="T") return true;
+    return false;
+}
+
+bool pFlow::isNo(const word & str)
+{
+    word s = toUpper(str);
+
+    if( s == "NO" || s == "N" || "FALSE" || s == "OFF" || s == "F") return true;
+    return false;
+}
+
+
+pFlow::word pFlow::angleBracketsNames(const word& w1, const word& w2)
+{
+    return w1+"<"+w2+">";
+}
+
+pFlow::word pFlow::angleBracketsNames2
+(
+    const word& base,
+    const word& w1,
+    const word& w2
+)
+{
+    return base+"<"+w1+","+w2+">";
+}
+
+pFlow::word pFlow::angleBracketsNames3(const word& base, const word& w1, const word& w2, const word& w3)
+{
+    return base+"<"+w1+","+w2+","+w3+">";
+}
+
+pFlow::word pFlow::groupNames(const word& bw, const word& tw, char sep)
+{
+    return bw + sep + tw;
+}
+
+pFlow::word pFlow::baseName(const word& w, char sep)
+{
+    if( auto pos = w.find_last_of(sep); pos != word::npos)
+    {
+        return w.substr(0,pos);
+    }
+    else
+    {
+        return w;
+    }
+}
+
+pFlow::word pFlow::tailName(const word& w, char sep)
+{
+    if( auto pos = w.find_last_of(sep); pos != word::npos)
+    {
+        return w.substr(pos+1);
+    }
+    else
+    {
+        return nullWord;
+    }   
+}
+
+bool pFlow::validWord(char c)
+{
+    return
+    (
+        !isspace(c)
+     && c != '"'   // string quote
+     && c != '\''  // string quote
+     //&& c != '/'   // path separator
+     && c != ';'   // end statement
+     && c != '{'   // beg subdict
+     && c != '}'   // end subdict
+    );
+}
+
+bool pFlow::validWord(const word& w)
+{
+    for(auto wi:w)
+    {   
+        char c = wi;
+        if ( !validWord(c) ) return false;
+    }
+    return true;
+}
+
+
+bool pFlow::readLabel( const word& w, label & val)
+{
+    try{
+        val = std::stoull(w);
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readLabel( const char* buf, label & val)
+{
+    word w(buf);
+    return readLabel(w, val);
+}
+
+bool pFlow::readUint32( const word& w, uint32 & val)
+{
+    try{
+        val = std::stoul(w);
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readUint32( const char* buf, uint32 & val)
+{
+    word w(buf);
+    return readUint32(w, val);
+}
+
+bool pFlow::readInt64( const word& w, int64 & val)
+{
+    try{
+        val = std::stoll(w);
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readInt64( const char* buf, int64 & val)
+{
+    word w(buf);
+    return readInt64(w, val);
+}
+
+bool pFlow::readInt32( const word& w, int32 & val)
+{
+    try{
+        val = std::stoi(w);
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readInt32( const char* buf, int32 & val)
+{
+    word w(buf);
+    return readInt32(w, val);
+}
+
+
+bool pFlow::readInt16( const word& w, int16 & val)
+{
+    try{
+        val = static_cast<int16>(std::stoi(w));
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readInt16( const char* buf, int16 & val)
+{
+    word w(buf);
+    return readInt16(w, val);
+}
+
+bool pFlow::readInt8( const word& w, int8 & val)
+{
+    try{
+        val = static_cast<int8>(std::stoi(w));
+    }
+    catch (...)
+    {
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readInt8( const char* buf, int8 & val)
+{
+    word w(buf);
+    return readInt8(w, val);
+}
+
+bool pFlow::readReal( const word& w, real & val)
+{
+    try{
+        val = std::stod(w);
+        
+    }
+    catch (...){
+        return false;
+    }
+    return true;
+}
+
+bool pFlow::readReal( const char* buf, real & val )
+{
+    word w(buf);
+    return readReal(w, val);
+}
+
+
+bool pFlow::readBoolian_Str( const word& w, bool & val)
+{
+    if( bool t = isYes(w); t )
+    {
+        val = true;
+        return true;    
+    } 
+    if( bool f = isNo(w); f )
+    {
+        val = false;
+        return true;
+    }
+    return false;
+}
+
+bool pFlow::readBoolian_Str( const char* buf, bool & val)
+{
+    word w(buf);
+    return readBoolian_Str(w, val);
+}
diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.H b/src/phasicFlow/types/basicTypes/bTypesFunctions.H
new file mode 100755
index 00000000..748b53f0
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.H
@@ -0,0 +1,236 @@
+/*------------------------------- 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 __bTypesFunctions_H__
+#define __bTypesFunctions_H__
+
+#include "dFlowMacros.H"
+#include "numericConstants.H"
+#include "builtinTypes.H"
+#include "math.H"
+
+
+// helper functions and constants for basic types 
+namespace pFlow
+{
+
+
+inline const real  zero = 0.0;
+inline const real  one  = 1.0;
+
+inline const int32  zero32 = 0;
+inline const int32  one32  = 1;
+
+inline const word nullWord;
+
+inline const word whiteSpace(" \t\n\v\f\r");
+
+
+int32 countChar (const word& s, const char c);
+
+int32 countChar(const char* s, const char c);
+
+word toUpper(const word & inStr);
+
+bool isYes(const word & str);
+
+bool isNo(const word & str);
+
+word real2Fixed(const real & v, int32 numPrecision = 6);
+
+word real2Word(const real & v, int32 numPrecision = 6);
+
+word int322Word(const int32 & v);
+
+ 
+word angleBracketsNames(const word& w1, const word& w2);
+
+ 
+word angleBracketsNames2(const word& base, const word& w1, const word& w2);
+
+word angleBracketsNames3(const word& base, const word& w1, const word& w2, const word& w3);
+
+word groupNames(const word& bw, const word& tw, char sep = '.');
+
+word baseName(const word& w, char sep = '.');
+
+word tailName(const word& w, char sep = '.');
+
+
+// is the character valid for a word name  
+bool validWord(char c);
+
+
+bool validWord(const word& w);
+
+
+
+bool readLabel( const word& w, label & val);
+
+
+bool readLabel( const char* buf, label & val);
+
+
+bool readUint32( const word& w, uint32 & val);
+
+
+bool readUint32( const char* buf, uint32 & val);
+
+
+bool readInt64( const word& w, int64 & val);
+
+
+bool readInt64( const char* buf, int64 & val);
+
+
+bool readInt32( const word& w, int32 & val);
+
+
+bool readInt32( const char* buf, int32 & val);
+
+
+bool readInt16( const word& w, int16 & val);
+
+
+bool readInt16( const char* buf, int16 & val);
+
+
+bool readInt8( const word& w, int8 & val);
+
+
+bool readInt8( const char* buf, int8 & val);
+
+
+bool readReal( const word& w, real & val);
+
+
+bool readReal( const char* buf, real & val );
+
+
+bool readBoolian_Str( const word& w, bool & val);
+
+
+bool readBoolian_Str( const char* buf, bool & val);
+
+inline
+bool readValue(const word& w, real& val)
+{
+    return readReal(w,val);
+}
+
+inline
+bool readValue(const word& w, label& val)
+{
+    return readLabel(w,val);
+}
+
+inline
+bool readValue(const word& w, uint32& val)
+{
+    return readUint32(w,val);
+}
+
+inline
+bool readValue(const word& w, int64& val)
+{
+    return readInt64(w,val);
+}
+
+inline
+bool readValue(const word& w, int32& val)
+{
+    return readInt32(w,val);
+}
+
+inline
+bool readValue(const word& w, int16& val)
+{
+    return readInt16(w,val);
+}
+
+inline
+bool readValue(const word& w, int8& val)
+{
+    return readInt8(w,val);
+}
+
+inline
+bool readValue(const word& w, bool& val)
+{
+    return readBoolian_Str(w,val);
+}
+
+
+INLINE_FUNCTION_HD bool equal(const real& s1, const real& s2)
+{
+    return abs(s1 - s2) <= smallValue;
+}
+
+INLINE_FUNCTION_HD bool equal(const int64& s1, const int64& s2)
+{
+    return s1 == s2;
+}
+
+INLINE_FUNCTION_HD bool equal(const int32& s1, const int32& s2)
+{
+	return s1 == s2;
+}
+
+INLINE_FUNCTION_HD bool equal(const int16& s1, const int16& s2)
+{
+    return s1 == s2;
+}
+
+INLINE_FUNCTION_HD bool equal(const int8& s1, const int8& s2)
+{
+    return s1 == s2;
+}
+
+INLINE_FUNCTION_HD bool equal(const uint32& s1, const uint32& s2)
+{
+    return s1 == s2;
+}
+
+INLINE_FUNCTION_HD bool equal(const label& s1, const label& s2)
+{
+    return s1 == s2;
+}
+
+// host only
+INLINE_FUNCTION bool equal(const word& s1, const word& s2)
+{
+	return s1==s2;
+}
+
+INLINE_FUNCTION_HD real degree2Radian(const real &theta)
+{
+    return theta / 180.0 * Pi;
+}
+
+INLINE_FUNCTION_HD real radian2Degree(const real &phi)
+{
+    return phi / Pi * 180.0;
+}
+
+
+} // end of pFlow
+
+
+#endif //__bTypesFunctions_H__
diff --git a/src/phasicFlow/types/basicTypes/builtinTypes.H b/src/phasicFlow/types/basicTypes/builtinTypes.H
new file mode 100755
index 00000000..19d38e60
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/builtinTypes.H
@@ -0,0 +1,57 @@
+/*------------------------------- 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 __builtinTypes_H__
+#define __builtinTypes_H__
+
+#include <string>
+
+
+namespace pFlow
+{
+
+#define useDouble 1
+
+// scalars
+#if useDouble
+	using real 	= double;
+#else
+	using real 	= float;
+#endif
+
+using int8 		= 	signed char;
+
+using int16 	=	short int;
+
+using int32 	=	int;
+
+using int64 	= 	long long int;
+
+using uint32 	= 	unsigned int;
+
+using label    	=  	std::size_t;
+
+using word 		= 	std:string;
+
+
+} // end of pFlow
+
+
+#endif
diff --git a/src/phasicFlow/types/basicTypes/math.H b/src/phasicFlow/types/basicTypes/math.H
new file mode 100755
index 00000000..875e5a9c
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/math.H
@@ -0,0 +1,410 @@
+/*------------------------------- 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 __math_H__
+#define __math_H__
+
+#include <cmath>
+
+#ifdef __CUDACC__
+#include "math.h"
+#endif
+
+#include "builtinTypes.H"
+
+//* * * * * * * * * * *  List of functinos * * * * * * * * //
+// abs, mod, exp, log, log10, pow, sqrt, cbrt
+// sin, cos, tan, asin, acos, atan, atan2
+// sinh, cosh, tanh, asinh, acosh, atanh 
+// min, max
+//* * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace dFlow
+{
+
+
+INLINE_FUNCTION_HD real abs(real x)
+{
+#ifdef __CUDACC__
+	return ::fabs(x);
+#else
+	return std::fabs(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD int64 abs(int64 x)
+{
+#ifdef __CUDACC__
+	return ::abs(x);
+#else
+	return std::abs(x);
+#endif
+}
+
+INLINE_FUNCTION_HD int32 abs(int32 x)
+{
+#ifdef __CUDACC__
+	return ::abs(x);
+#else
+	return std::abs(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real mod(real x, real y)
+{
+#ifdef __CUDACC__
+	return ::fmod(x, y);
+#else
+	return std::fmod(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD int64 mod(int64 x, int64 y)
+{
+	return x%y;
+}
+
+INLINE_FUNCTION_HD int32 mod(int32 x, int32 y)
+{
+	return x%y;
+}
+
+INLINE_FUNCTION_HD int64 mod(label x, label y)
+{
+	return x%y;
+}
+
+INLINE_FUNCTION_HD int32 mod(uint32 x, uint32 y)
+{
+	return x%y;
+}
+
+INLINE_FUNCTION_HD real remainder(real x, real y)
+{
+#ifdef __CUDACC__
+	return ::remainder(x,y);
+#else
+	return std::remainder(x,y);
+#endif
+}
+
+INLINE_FUNCTION_HD real exp(real x)
+{
+#ifdef __CUDACC__
+	return ::exp(x);
+#else
+	return std::exp(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real log(real x)
+{
+#ifdef __CUDACC__
+	return ::log(x);
+#else
+	return std::log(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real log10(real x)
+{
+#ifdef __CUDACC__
+	return ::log10(x);
+#else
+	return std::log10(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real pow(real x, real y)
+{
+#ifdef __CUDACC__
+	return ::pow(x, y);
+#else
+	return std::pow(x, y);
+#endif
+
+}
+
+
+INLINE_FUNCTION_HD real sqrt(real x)
+{
+#ifdef __CUDACC__
+	return ::sqrt(x);
+#else
+	return std::sqrt(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real cbrt (real x)
+{
+#ifdef __CUDACC__
+	return ::cbrt (x);
+#else
+	return std::cbrt (x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real sin(real x)
+{
+#ifdef __CUDACC__
+	return ::sin(x);
+#else
+	return std::sin(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real cos(real x)
+{
+#ifdef __CUDACC__
+	return ::cos(x);
+#else
+	return std::cos(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real tan(real x)
+{
+#ifdef __CUDACC__
+	return ::tan(x);
+#else
+	return std::tan(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real asin(real x)
+{
+#ifdef __CUDACC__
+	return ::asin(x);
+#else
+	return std::asin(x);
+#endif
+	
+}
+
+
+INLINE_FUNCTION_HD real acos(real x)
+{
+#ifdef __CUDACC__
+	return ::acos(x);
+#else
+	return std::acos(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real atan(real x)
+{
+#ifdef __CUDACC__
+	return ::atan(x);
+#else
+	return std::atan(x);
+#endif
+
+}
+
+
+INLINE_FUNCTION_HD real atan2(real y, real x)
+{
+#ifdef __CUDACC__
+	return ::atan2(y, x);
+#else
+	return std::atan2(y, x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real sinh(real x)
+{
+#ifdef __CUDACC__
+	return ::sinh(x);
+#else
+	return std::sinh(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real cosh(real x)
+{
+#ifdef __CUDACC__
+	return ::cosh(x);
+#else
+	return std::cosh(x);
+#endif
+}
+
+
+INLINE_FUNCTION_HD real tanh(real x)
+{
+#ifdef __CUDACC__
+	return ::tanh(x);
+#else
+	return std::tanh(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real asinh(real x)
+{
+#ifdef __CUDACC__
+	return ::asinh(x);
+#else
+	return std::asinh(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real acosh(real x)
+{
+#ifdef __CUDACC__
+	return ::acosh(x);
+#else
+	return std::acosh(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real atanh(real x)
+{
+#ifdef __CUDACC__
+	return ::atanh(x);
+#else
+	return std::atanh(x);
+#endif
+}
+
+INLINE_FUNCTION_HD real min(real x, real y)
+{
+#ifdef __CUDACC__
+	return ::fmin(x,y);
+#else
+	return std::fmin(x,y);
+#endif
+}
+
+INLINE_FUNCTION_HD int64 min(int32 x, int32 y)
+{
+#ifdef __CUDACC__
+	return ::min(x, y);
+#else
+	return std::min(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD int64 min(int64 x, int64 y)
+{
+#ifdef __CUDACC__
+	return ::min(x, y);
+#else
+	return std::min(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD label min(label x, label y)
+{
+#ifdef __CUDACC__
+	return ::min(x, y);
+#else
+	return std::min(x, y);
+#endif
+
+}
+
+INLINE_FUNCTION_HD uint32 min(uint32 x, uint32 y)
+{
+#ifdef __CUDACC__
+	return ::min(x, y);
+#else
+	return std::min(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD uint32 min(uint16 x, uint16 y)
+{
+#ifdef __CUDACC__
+	return ::min(x, y);
+#else
+	return std::min(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD real max(real x, real y)
+{
+#ifdef __CUDACC__
+	return ::fmax(x,y);
+#else
+	return std::fmax(x,y);
+#endif
+}
+
+INLINE_FUNCTION_HD int64 max(int64 x, int64 y)
+{
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
+	return std::max(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD int32 max(int32 x, int32 y)
+{
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
+	return std::max(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD label max(label x, label y)
+{
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
+	return std::max(x, y);
+#endif
+}
+
+
+INLINE_FUNCTION_HD uint32 max(uint32 x, uint32 y)
+{
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
+	return std::max(x, y);
+#endif
+}
+
+INLINE_FUNCTION_HD uint32 max(uint16 x, uint16 y)
+{
+#ifdef __CUDACC__
+	return ::max(x, y);
+#else
+	return std::max(x, y);
+#endif
+
+}
+
+} // pFlow
+
+
+
+#endif // __math_H__
diff --git a/src/phasicFlow/types/basicTypes/numericConstants.H b/src/phasicFlow/types/basicTypes/numericConstants.H
new file mode 100755
index 00000000..2fba1538
--- /dev/null
+++ b/src/phasicFlow/types/basicTypes/numericConstants.H
@@ -0,0 +1,72 @@
+/*------------------------------- 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 __numericConstants_H__
+#define __numericConstants_H__
+
+
+#include <limits>
+#include "builtinTypes.H"
+
+
+namespace pFlow
+{
+
+	const inline real Pi = real(3.1415926535897932384626433832);
+	const inline real smallValue = 1.0e-15;
+	const inline real verySmallValue = 1.0e-30;
+	const inline real largeValue = 1.0e15;
+	const inline real veryLargeValue = 1.0e30;
+
+	// - largest negative value 
+	template<typename T>
+	constexpr inline T largestNegative()
+	{
+		return std::numeric_limits<T>::lowest();
+	}
+
+	template<typename T>
+	constexpr inline T epsilonValue()
+	{
+		return std::numeric_limits<T>::min();
+	}
+
+	// largest positive value 
+	template<typename T>
+	constexpr inline T largestPositive()
+	{
+		return std::numeric_limits<T>::max();
+	}
+
+	const inline int32 largestNegInt32  		= largestNegative<int32>();
+	const inline int32 largestPosInt32  		= largestPositive<int32>();
+
+	const inline int64 largestNegInt64  		= largestNegative<int64>();
+	const inline int64 largestPosInt64  		= largestPositive<int64>();
+		
+	const inline real largestNegREAL 	= largestNegative<real>();
+	const inline real largestPosREAL 	= largestPositive<real>();
+	const inline real epsilonREAL 		= epsilonValue<real>();
+	
+	
+
+} // end of pFlow
+
+#endif //__numericConstants_H__
diff --git a/src/phasicFlow/types/quadruple/quadruple.H b/src/phasicFlow/types/quadruple/quadruple.H
new file mode 100644
index 00000000..7263d76a
--- /dev/null
+++ b/src/phasicFlow/types/quadruple/quadruple.H
@@ -0,0 +1,240 @@
+/*------------------------------- 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 __quadruple_H__
+#define __quadruple_H__
+
+
+#include "uniquePtr.H"
+#include "triple.H"
+#include "iOstream.H"
+#include "iIstream.H"
+#include "token.H"
+#include "error.H"
+
+
+namespace pFlow
+{
+
+
+template<typename T> class quadruple;
+class iIstream;
+
+#include "quadrupleFwd.H"
+
+// you can see it as a sequence of four elements (w,x,y,z) or an (scalar and vector)
+template<typename T>
+class quadruple
+{
+public:
+
+	T 		  s_;
+	triple<T> v_;
+public:
+	
+	//// constructors
+	INLINE_FUNCTION_HD
+	quadruple()
+	{}
+
+	//// Constructors
+
+	INLINE_FUNCTION_HD
+	quadruple(const T &w, const T &x, const T &y, const T &z)
+	:
+		s_(w),
+		v_(x, y, z)
+	{}
+
+	INLINE_FUNCTION_HD
+	quadruple(const T& s, const triple<T> v)
+	:
+		s_(s),
+		v_(v)
+	{}
+
+	INLINE_FUNCTION_HD
+	quadruple(const T &val)
+	:
+		s_(val),
+		v_(val)
+	{}
+
+	// type conversion trough assignment  
+	template <typename T2>
+	INLINE_FUNCTION_HD
+	quadruple<T> & operator = (const quadruple<T2> & rhs)
+	{
+		this->v_ = rhs.v_;
+		this->s_ = static_cast<T>(rhs.s_);
+		return *this;
+	}
+
+	// type casting through copy constructor 
+	template<typename T2>
+	INLINE_FUNCTION_HD
+	quadruple(const quadruple<T2> &src):
+	s_(static_cast<T>(src.s_)),
+	v_(src.v_)
+	{}
+
+	// copy construct
+	INLINE_FUNCTION_HD
+	quadruple(const quadruple<T>& src) = default;
+
+	// move construct
+	INLINE_FUNCTION_HD
+	quadruple(quadruple<T>&& src) = default;
+
+	// copy assignment
+	INLINE_FUNCTION_HD
+	quadruple<T>& operator=(const quadruple<T>& src) = default;
+
+	// move assignment
+	INLINE_FUNCTION_HD
+	quadruple<T>& operator=(quadruple<T>&& src) = default;
+
+	// clone 
+	INLINE_FUNCTION_H
+	uniquePtr<quadruple<T>> clone()const
+	{
+		return makeUnique<quadruple<T>>(*this);
+	}
+	
+	INLINE_FUNCTION_H
+	quadruple<T>* clonePtr()const
+	{
+		return new quadruple<T>(*this);
+	}
+
+	// Access
+	INLINE_FUNCTION_HD
+	T & w(){ return s_; }
+	INLINE_FUNCTION_HD
+	const T & w()const { return s_; }
+
+	INLINE_FUNCTION_HD
+	T & x(){ return v_.x(); }
+	INLINE_FUNCTION_HD
+	const T & x()const { return v_.x(); }
+	
+	INLINE_FUNCTION_HD
+	T & y(){ return v_.y(); }
+	INLINE_FUNCTION_HD
+	const T & y()const { return v_.y(); }
+
+	INLINE_FUNCTION_HD
+	T & z(){ return v_.z(); }
+	INLINE_FUNCTION_HD
+	const T & z()const { return v_.z(); }
+	
+	INLINE_FUNCTION_HD
+	T & s(){ return s_; }
+	INLINE_FUNCTION_HD
+	const T & s()const { return s_; }	
+	
+	INLINE_FUNCTION_HD
+	triple<T> v() {return v_;}
+	INLINE_FUNCTION_HD
+	const triple<T> v() const {return v_;}
+
+	
+	// methods
+	friend FUNCTION_HD T dot <T> (const quadruple<T> & oprnd1, const quadruple<T> & oprnd2);
+	
+	INLINE_FUNCTION_HD T length() const; 
+
+	INLINE_FUNCTION_HD void normalize(); 
+
+
+	//// operators 
+	
+	// + operator
+	friend FUNCTION_HD quadruple<T> operator+ <T> (const quadruple<T> & oprnd1, const quadruple<T> & oprnd2);
+
+	friend FUNCTION_HD quadruple<T> operator+ <T> (const quadruple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator+ <T> (const T & oprnd1, const quadruple<T> & oprnd2);	
+	
+	// - operator 
+	friend FUNCTION_HD quadruple<T> operator - <T> (const quadruple<T> & oprnd1, const quadruple<T> & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator - <T> (const quadruple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator - <T> (const T & oprnd1, const quadruple<T> & oprnd2);
+	
+	// * operators
+	friend FUNCTION_HD quadruple<T> operator * <T> (const quadruple<T> & oprnd1, const quadruple<T> & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator * <T> (const quadruple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator * <T> (const T & oprnd1, const quadruple<T> & oprnd2);
+
+	// / operators
+	friend FUNCTION_HD quadruple<T> operator / <T> (const quadruple<T> & oprnd1, const quadruple<T> & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator / <T> (const quadruple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD quadruple<T> operator / <T> (const T & oprnd1, const quadruple<T> & oprnd2);
+
+
+	INLINE_FUNCTION_HD
+	void operator+= (const quadruple & oprnd2);
+
+	INLINE_FUNCTION_HD
+	void operator-= (const quadruple & oprnd2);
+
+	INLINE_FUNCTION_HD
+	void operator*= (const quadruple & oprnd2);
+
+	INLINE_FUNCTION_HD
+	void operator/= (const quadruple & oprnd2);
+	
+	
+	// unary negate operator 
+	INLINE_FUNCTION_HD
+	quadruple operator- ()const;
+	
+	// unary plus operator 
+	INLINE_FUNCTION_HD
+	quadruple operator+ ()const;
+	
+
+	friend FUNCTION_HD bool operator == <T> (const quadruple<T> &opr1, const quadruple<T> &opr2);
+	
+	// << operator 	
+	friend FUNCTION_H iOstream& operator<< <T> (iOstream& str, const quadruple<T> & ov);
+	
+	// >> operator 
+	friend FUNCTION_H iIstream& operator >> <T> (iIstream & str, quadruple<T> & iv);
+
+	// same as >> operator, but faster, good for mass read
+	friend FUNCTION_H void readIstream <T>( iIstream& str, quadruple<T> &iv);
+
+
+};
+
+} // pFlow
+
+
+#include "quadrupleI.H"
+//#include "quadrupleMath.H"
+
+#endif
diff --git a/src/phasicFlow/types/quadruple/quadrupleFwd.H b/src/phasicFlow/types/quadruple/quadrupleFwd.H
new file mode 100644
index 00000000..7ba0a085
--- /dev/null
+++ b/src/phasicFlow/types/quadruple/quadrupleFwd.H
@@ -0,0 +1,141 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename T>
+INLINE_FUNCTION_HD T dot
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+);
+
+
+template<typename T> 
+INLINE_FUNCTION_HD quadruple<T> operator + 
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator+
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+);
+
+
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator+
+(
+	const T & oprnd2,
+	const quadruple<T> & oprnd1
+); 
+
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator-
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator-
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator-
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator*
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator*
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator*
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator/
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator/
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD quadruple<T> operator/
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator ==
+(
+	const quadruple<T> &opr1,
+	const quadruple<T> &opr2
+);
+
+template<typename T> 
+INLINE_FUNCTION_H iOstream& operator<<
+(
+	iOstream& str,
+	const quadruple<T> & ov
+);
+
+
+template<typename T> 
+INLINE_FUNCTION_H iIstream& operator>>
+(
+	iIstream& str,
+	quadruple<T> & iv
+);
+
+template<typename T>
+INLINE_FUNCTION_H void readIstream 
+(
+	iIstream& str,
+	quadruple<T> & iv
+);
\ No newline at end of file
diff --git a/src/phasicFlow/types/quadruple/quadrupleI.H b/src/phasicFlow/types/quadruple/quadrupleI.H
new file mode 100644
index 00000000..7c45094f
--- /dev/null
+++ b/src/phasicFlow/types/quadruple/quadrupleI.H
@@ -0,0 +1,351 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename T> 
+INLINE_FUNCTION_HD  T pFlow::dot
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return  oprnd1.s_ * oprnd2.s_ + 
+ 			dot(oprnd1.v(), oprnd2.v()) ;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD T pFlow::quadruple<T>::length
+(
+)const
+{
+	return sqrt(dot(*this,*this));
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::quadruple<T>::normalize
+(
+)
+{
+	T l = length();
+	if( static_cast<real>(l) > smallValue ) 
+	{
+		*this /= l;
+	}
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator+ 
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ + oprnd2.s_,
+		oprnd1.v_ + oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator+
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ + oprnd2,
+		oprnd1.v_ + oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator+
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1 + oprnd2.s_,
+		oprnd1 + oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator- 
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ - oprnd2.s_,
+		oprnd1.v_ - oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator-
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ - oprnd2,
+		oprnd1.v_ - oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator-
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1 - oprnd2.s_,
+		oprnd1 - oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator*
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ * oprnd2.s_,
+		oprnd1.v_ * oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator*
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ * oprnd2,
+		oprnd1.v_ * oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator*
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1 * oprnd2.s_,
+		oprnd1 * oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator/
+(
+	const quadruple<T> & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ / oprnd2.s_,
+		oprnd1.v_ / oprnd2.v_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator/
+(
+	const quadruple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1.s_ / oprnd2,
+		oprnd1.v_ / oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::operator/
+(
+	const T & oprnd1,
+	const quadruple<T> & oprnd2
+)
+{
+	return quadruple<T>(
+		oprnd1 / oprnd2.s_,
+		oprnd1 / oprnd2.v_
+		);
+}
+ 
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::quadruple<T>::operator+=
+(
+	const quadruple<T> & oprnd2
+)
+{	
+	this->s_ += oprnd2.s_;
+	this->v_ += oprnd2.v_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::quadruple<T>::operator-=
+(
+	const quadruple<T> & oprnd2
+)
+{	
+	this->s_ -= oprnd2.s_;
+	this->v_ -= oprnd2.v_;
+	
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::quadruple<T>::operator*=
+(
+	const quadruple<T> & oprnd2
+)
+{	
+	this->s_ *= oprnd2.s_;
+	this->v_ *= oprnd2.v_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::quadruple<T>::operator/=
+(
+	const quadruple<T> & oprnd2
+)
+{	
+	this->s_ /= oprnd2.s_;
+	this->v_ /= oprnd2.v_;
+	
+}
+
+template<typename T>
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::quadruple<T>::operator-
+(
+) const
+{
+	return quadruple<T>(-this->s_, -this->v_);
+}
+	
+template<typename T>
+INLINE_FUNCTION_HD pFlow::quadruple<T> pFlow::quadruple<T>::operator+
+(
+) const
+{
+	return *this;
+}
+
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator == 
+(
+	const quadruple<T> &opr1,
+	const quadruple<T> &opr2
+){
+	return equal(opr1.s_, opr2.s_) && equal(opr1.v_, opr2.v_);
+};
+
+template<typename T>
+INLINE_FUNCTION_H pFlow::iOstream& pFlow::operator <<
+(
+	iOstream & str,
+	const quadruple<T> & ov
+)
+{
+
+	str << token::BEGIN_LIST << ov.w()
+	    << token::SPACE << ov.x()
+	    << token::SPACE << ov.y()
+	    << token::SPACE << ov.z()
+		<< token::END_LIST;	    
+
+	str.check(FUNCTION_NAME);
+
+	return str;
+}
+
+template<typename T>
+INLINE_FUNCTION_H pFlow::iIstream& pFlow::operator >>
+(
+	iIstream & str,
+	quadruple<T> & iv
+)
+{
+
+	str.readBegin("quadruple<T>");
+
+	str >> iv.w();
+    str >> iv.x();
+    str >> iv.y();
+    str >> iv.z();
+
+    str.readEnd("quadruple<T>");
+
+    str.check(FUNCTION_NAME);
+
+	return str;
+}
+
+template<typename T>
+INLINE_FUNCTION_H void pFlow::readIstream
+(
+	iIstream & str,
+	quadruple<T> & iv
+)
+{
+
+	str.readBegin("quadruple<T>");
+	T val;
+
+    readIstream(str, val);
+    iv.w() = val;
+
+    readIstream(str, val);
+    iv.x() = val;
+    
+    readIstream(str, val);
+    iv.y() = val;
+
+    readIstream(str, val);
+    iv.z() = val;
+    
+    str.readEnd("quadruple<T>");
+
+    str.check(FUNCTION_NAME);
+
+}
diff --git a/src/phasicFlow/types/quadruple/quadrupleMath.H b/src/phasicFlow/types/quadruple/quadrupleMath.H
new file mode 100644
index 00000000..68980765
--- /dev/null
+++ b/src/phasicFlow/types/quadruple/quadrupleMath.H
@@ -0,0 +1,91 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#define Q4Func(fnName) 									\
+template<typename T> 									\
+inline pFlow::quadruple<T> pFlow::fnName(const quadruple<T>& q) 			\
+{																\
+	return quadruple<T>(fnName(q.s_), fnName(q.v_));	\
+}
+
+#define Q4Func2(fnName)  															\
+template<typename T>																\
+inline pFlow::quadruple<T> pFlow::fnName(const quadruple<T>& arg1, const quadruple<T>& arg2)	\
+{																					\
+	return quadruple<T>(fnName(arg1.s_, arg2.s_), fnName(arg1.v_,arg2.v_)); \
+}
+
+//* * * * * * * * * * *  List of functinos * * * * * * * * //
+// abs, mod, exp, log, log10, pow, sqrt, cbrt
+// sin, cos, tan, asin, acos, atan, atan2
+// sinh, cosh, tanh, asinh, acosh, atanh 
+// min, max
+//* * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+Q4Func(abs);
+Q4Func2(mod);
+Q4Func(exp);
+Q4Func(log);
+Q4Func(log10);
+Q4Func2(pow);
+Q4Func(sqrt);
+Q4Func(cbrt);
+Q4Func(sin);
+Q4Func(cos);
+Q4Func(tan);
+Q4Func(asin);
+Q4Func(acos);
+Q4Func(atan);
+Q4Func2(atan2);
+Q4Func(sinh);
+Q4Func(cosh);
+Q4Func(tanh);
+Q4Func(asinh);
+Q4Func(acosh);
+Q4Func(atanh);
+Q4Func2(min);
+Q4Func2(max);
+
+
+template<typename T>
+inline pFlow::quadruple<T> pFlow::pow(const quadruple<T>& q4, T e)
+{
+	return quadruple<T>( pow(q4.s_, e), pow(q4.v_,e));
+}
+
+// return the min of 3 elements x, y, z
+template<typename T>
+inline T pFlow::min(const quadruple<T>& q4)
+{
+	return min( min(q4.v_), q4.s_);
+}
+
+// return the max of 3 elements x, y, z
+template<typename T>
+inline T pFlow::max(const quadruple<T>& q4)
+{
+	return max( max(q4.v_), q4.s_);
+}
+
+
+
+#undef Q4Func
+#undef Q4Func2
\ No newline at end of file
diff --git a/src/phasicFlow/types/triple/triple.H b/src/phasicFlow/types/triple/triple.H
new file mode 100644
index 00000000..b2029080
--- /dev/null
+++ b/src/phasicFlow/types/triple/triple.H
@@ -0,0 +1,236 @@
+/*------------------------------- 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 __triple_H__
+#define __triple_H__
+
+#include "pFlowMacros.H"
+#include "numericConstants.H"
+#include "uniquePtr.H"
+#include "iOstream.H"
+#include "iIstream.H"
+#include "error.H"
+
+
+
+namespace pFlow
+{
+
+
+template<typename T> class triple;
+
+
+#include "tripleFwd.H"
+
+
+// for 3D vectors
+// it should be used only for numeric types, real, unit 
+template <typename T>
+struct triple
+{
+	// data members
+	T x_;
+	T y_;
+	T z_;
+
+	// initilizes to zero
+	INLINE_FUNCTION_HD triple():
+	x_(0),
+	y_(0),
+	z_(0)
+	{}
+
+	// Constructors
+	INLINE_FUNCTION_HD triple(const T &x, const T &y, const T &z):
+	x_(x),
+	y_(y),
+	z_(z)
+	{}
+
+	INLINE_FUNCTION_HD triple(const T &v):
+	triple(v, v, v)
+	{}
+
+	// type conversion trough assignment  
+	template <typename T2>
+	INLINE_FUNCTION_HD triple<T> & operator = (const triple<T2> & rhs)
+	{
+		this->x_ = static_cast<T>(rhs.x_);
+		this->y_ = static_cast<T>(rhs.y_);
+		this->z_ = static_cast<T>(rhs.z_);
+		return *this;
+	}
+
+	// type casting through copy constructor 
+	template<typename T2>
+	INLINE_FUNCTION_HD triple(const triple<T2> &src):
+	x_(static_cast<T>(src.x_)),
+	y_(static_cast<T>(src.y_)),
+	z_(static_cast<T>(src.z_))
+	{
+	}
+	
+	// copy construct
+	INLINE_FUNCTION_HD
+	triple(const triple<T>& src) = default;
+
+	// volatile copy construct
+	/*INLINE_FUNCTION_HD
+	triple(volatile triple<T>& src):
+		x_(src.x_),
+		y_(src.y_),
+		z_(src.z_)
+	{}*/
+
+	/*INLINE_FUNCTION_HD
+	triple& operator=(volatile triple<T>& src)	
+	{
+		x_ = src.x_;
+		y_ = src.y_;
+		z_ = src.z_;
+	}*/
+
+	// move construct
+	INLINE_FUNCTION_HD
+	triple(triple<T>&& src) = default;
+
+	// copy assignment
+	INLINE_FUNCTION_HD
+	triple<T>& operator=(const triple<T>& src) = default;
+
+	// move assignment
+	INLINE_FUNCTION_HD
+	triple<T>& operator=(triple<T>&& src) = default;
+
+	// clone 
+	INLINE_FUNCTION 
+	uniquePtr<triple<T>> clone() const
+	{
+		return makeUnique<triple<T>>(*this);
+	}
+
+	INLINE_FUNCTION
+	triple<T>* clonePtr()const
+	{
+		return new triple<T>(*this);
+	}
+
+	//// member methods
+
+	// access
+	INLINE_FUNCTION_HD T & x(){ return x_; }
+	INLINE_FUNCTION_HD const T & x()const { return x_; }
+	
+	INLINE_FUNCTION_HD T & y(){ return y_; }
+	INLINE_FUNCTION_HD const T & y()const { return y_; }
+
+	INLINE_FUNCTION_HD T & z(){ return z_; }
+	INLINE_FUNCTION_HD const T & z()const { return z_; }
+
+	// methods
+	friend FUNCTION_HD T dot <T> (const triple<T> & oprnd1, const triple<T> & oprnd2);
+	
+	friend FUNCTION_HD triple<T> cross <T>(const triple<T> & v1, const triple<T> & v2);
+	
+	INLINE_FUNCTION_HD T length() const; 
+
+	INLINE_FUNCTION_HD void normalize(); 
+
+
+
+	//// operators 
+	
+	// + operator
+	friend FUNCTION_HD triple<T> operator+ <T> (const triple<T> & oprnd1, const triple<T> & oprnd2);
+
+	friend FUNCTION_HD triple<T> operator+ <T> (const triple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator+ <T> (const T & oprnd1, const triple<T> & oprnd2);	
+	
+	// - operator 
+	friend FUNCTION_HD triple<T> operator - <T> (const triple<T> & oprnd1, const triple<T> & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator - <T> (const triple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator - <T> (const T & oprnd1, const triple<T> & oprnd2);
+	
+	// * operators
+	friend FUNCTION_HD triple<T> operator * <T> (const triple<T> & oprnd1, const triple<T> & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator * <T> (const triple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator * <T> (const T & oprnd1, const triple<T> & oprnd2);
+
+	// / operators
+	friend FUNCTION_HD triple<T> operator / <T> (const triple<T> & oprnd1, const triple<T> & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator / <T> (const triple<T> & oprnd1, const T & oprnd2);
+	
+	friend FUNCTION_HD triple<T> operator / <T> (const T & oprnd1, const triple<T> & oprnd2);
+
+
+	INLINE_FUNCTION_HD void operator+= (const triple & oprnd2);
+
+	INLINE_FUNCTION_HD void operator-= (const triple & oprnd2);
+
+	INLINE_FUNCTION_HD void operator*= (const triple & oprnd2);
+
+	INLINE_FUNCTION_HD void operator/= (const triple & oprnd2);
+	
+	
+	// unary negate operator 
+	INLINE_FUNCTION_HD triple operator- ()const;
+	
+	// unary plus operator 
+	INLINE_FUNCTION_HD triple operator+ ()const;
+	
+
+	friend FUNCTION_HD bool operator == <T> (const triple<T> &opr1, const triple<T> &opr2);
+
+	friend FUNCTION_HD bool operator < <T> (const triple<T> &opr1, const triple<T> &opr2);
+
+	friend FUNCTION_HD bool operator > <T> (const triple<T> &opr1, const triple<T> &opr2);
+
+	friend FUNCTION_HD bool operator >= <T> (const triple<T> &opr1, const triple<T> &opr2);
+
+	friend FUNCTION_HD bool operator <= <T> (const triple<T> &opr1, const triple<T> &opr2);
+	
+	// << operator 	
+	friend iOstream& operator<< <T> (iOstream& str, const triple<T> & ov);
+	
+	// >> operator 
+	friend iIstream& operator >> <T> (iIstream & str, triple<T> & iv);
+
+	// same as >> operator, but faster, good for mass read
+	friend void readIstream <T>( iIstream& str, triple<T> &iv);
+
+};
+
+template<typename T>
+bool INLINE_FUNCTION_HD equal( const triple<T>& opr1, const triple<T>& opr2 );
+
+
+} // end of pFlow
+
+#include "tripleI.H"
+#include "tripleMath.H"
+
+
+#endif
diff --git a/src/phasicFlow/types/triple/tripleFwd.H b/src/phasicFlow/types/triple/tripleFwd.H
new file mode 100644
index 00000000..fe834702
--- /dev/null
+++ b/src/phasicFlow/types/triple/tripleFwd.H
@@ -0,0 +1,189 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename T>
+INLINE_FUNCTION_HD T dot
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> cross
+(
+	const triple<T> & v1,
+	const triple<T> & v2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD T length
+(
+	const triple<T> & v1
+);
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> normalize
+(
+	const triple<T>& v1
+);
+
+
+template<typename T> 
+INLINE_FUNCTION_HD triple<T> operator + 
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator+
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+);
+
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator+
+(
+	const T & oprnd2,
+	const triple<T> & oprnd1
+); 
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator-
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator-
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator-
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator*
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator*
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator*
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator/
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator/
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+);
+	
+template<typename T>
+INLINE_FUNCTION_HD triple<T> operator/
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator ==
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator >
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator <
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+);
+
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator <=
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+);
+
+template<typename T>
+INLINE_FUNCTION_HD bool operator >=
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+);
+
+template<typename T> 
+INLINE_FUNCTION iOstream& operator<<
+(
+	iOstream& str,
+	const triple<T> & ov
+);
+
+
+template<typename T> 
+INLINE_FUNCTION iIstream& operator>>
+(
+	iIstream& str,
+	triple<T> & iv
+);
+
+template<typename T>
+INLINE_FUNCTION void readIstream 
+(
+	iIstream& str,
+	triple<T> & iv
+);
\ No newline at end of file
diff --git a/src/phasicFlow/types/triple/tripleI.H b/src/phasicFlow/types/triple/tripleI.H
new file mode 100644
index 00000000..7d9ce5fe
--- /dev/null
+++ b/src/phasicFlow/types/triple/tripleI.H
@@ -0,0 +1,467 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename T> 
+INLINE_FUNCTION_HD  T pFlow::dot
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return  oprnd1.x_ * oprnd2.x_ + 
+ 			oprnd1.y_ * oprnd2.y_ + 
+ 			oprnd1.z_ * oprnd2.z_ ;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::cross
+(
+	const triple<T> & v1,
+	const triple<T> & v2
+)
+{
+	return triple<T>(
+			v1.y_*v2.z_ - v1.z_*v2.y_,
+			v1.z_*v2.x_ - v1.x_*v2.z_,
+			v1.x_*v2.y_ - v1.y_*v2.x_
+			);
+}
+
+template<typename T>
+INLINE_FUNCTION_HD T pFlow::length
+(
+	const triple<T> & v1
+)
+{
+	return v1.length();
+}
+
+template<typename T>
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::normalize(const triple<T>& v1)
+{
+	return v1/max(length(v1),verySmallValue);
+}
+
+template<typename T>
+INLINE_FUNCTION_HD T pFlow::triple<T>::length
+(
+)const
+{
+	return sqrt(dot(*this,*this));
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::triple<T>::normalize
+(
+)
+{
+	*this = *this/max(length(),verySmallValue);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator+ 
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ + oprnd2.x_,
+		oprnd1.y_ + oprnd2.y_,
+		oprnd1.z_ + oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator+
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ + oprnd2,
+		oprnd1.y_ + oprnd2,
+		oprnd1.z_ + oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator+
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1 + oprnd2.x_,
+		oprnd1 + oprnd2.y_,
+		oprnd1 + oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator- 
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ - oprnd2.x_,
+		oprnd1.y_ - oprnd2.y_,
+		oprnd1.z_ - oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator-
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ - oprnd2,
+		oprnd1.y_ - oprnd2,
+		oprnd1.z_ - oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator-
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1 - oprnd2.x_,
+		oprnd1 - oprnd2.y_,
+		oprnd1 - oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator*
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ * oprnd2.x_,
+		oprnd1.y_ * oprnd2.y_,
+		oprnd1.z_ * oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator*
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ * oprnd2,
+		oprnd1.y_ * oprnd2,
+		oprnd1.z_ * oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator*
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1 * oprnd2.x_,
+		oprnd1 * oprnd2.y_,
+		oprnd1 * oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator/
+(
+	const triple<T> & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ / oprnd2.x_,
+		oprnd1.y_ / oprnd2.y_,
+		oprnd1.z_ / oprnd2.z_
+		);
+}
+
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator/
+(
+	const triple<T> & oprnd1,
+	const T & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1.x_ / oprnd2,
+		oprnd1.y_ / oprnd2,
+		oprnd1.z_ / oprnd2
+		);
+}
+
+template<typename T> 
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::operator/
+(
+	const T & oprnd1,
+	const triple<T> & oprnd2
+)
+{
+	return triple<T>(
+		oprnd1 / oprnd2.x_,
+		oprnd1 / oprnd2.y_,
+		oprnd1 / oprnd2.z_
+		);
+}
+ 
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::triple<T>::operator+=
+(
+	const triple<T> & oprnd2
+)
+{	
+	this->x_ = this->x_ + oprnd2.x_;
+	this->y_ = this->y_ + oprnd2.y_;
+	this->z_ = this->z_ + oprnd2.z_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::triple<T>::operator-=
+(
+	const triple<T> & oprnd2
+)
+{	
+	this->x_ = this->x_ - oprnd2.x_;
+	this->y_ = this->y_ - oprnd2.y_;
+	this->z_ = this->z_ - oprnd2.z_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::triple<T>::operator*=
+(
+	const triple<T> & oprnd2
+)
+{	
+	this->x_ = this->x_ * oprnd2.x_;
+	this->y_ = this->y_ * oprnd2.y_;
+	this->z_ = this->z_ * oprnd2.z_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD void pFlow::triple<T>::operator/=
+(
+	const triple<T> & oprnd2
+)
+{	
+	this->x_ = this->x_ / oprnd2.x_;
+	this->y_ = this->y_ / oprnd2.y_;
+	this->z_ = this->z_ / oprnd2.z_;
+}
+
+template<typename T>
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::triple<T>::operator-
+(
+) const
+{
+	return triple<T>(-this->x_, -this->y_, -this->z_);
+}
+	
+template<typename T>
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::triple<T>::operator+
+(
+) const
+{
+	return *this;
+}
+
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator == 
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+){
+	return equal(opr1, opr2);
+};
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator <
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+)
+{
+	if( opr1.x_ < opr2.x_ && opr1.y_ < opr2.y_ && opr1.z_ < opr2.z_)
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+};
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator >
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+)
+{
+	if( opr1.x_ > opr2.x_ && opr1.y_ > opr2.y_ && opr1.z_ > opr2.z_)
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+};
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator <=
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+)
+{
+	if( opr1.x_ <= opr2.x_ && opr1.y_ <= opr2.y_ && opr1.z_ <= opr2.z_)
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+}
+
+
+template<typename T>
+INLINE_FUNCTION_HD bool pFlow::operator >=
+(
+	const triple<T> &opr1,
+	const triple<T> &opr2
+)
+{
+	if( opr1.x_ >= opr2.x_ && opr1.y_ >= opr2.y_ && opr1.z_ >= opr2.z_)
+	{
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+}
+
+
+template<typename T>
+INLINE_FUNCTION pFlow::iOstream& pFlow::operator <<
+(
+	iOstream & str,
+	const triple<T> & ov
+)
+{
+
+	str << token::BEGIN_LIST << ov.x_
+	    << token::SPACE << ov.y_
+	    << token::SPACE << ov.z_
+		<< token::END_LIST;	    
+
+	str.check(FUNCTION_NAME);
+
+	return str;
+}
+
+template<typename T>
+INLINE_FUNCTION pFlow::iIstream& pFlow::operator >>
+(
+	iIstream & str,
+	triple<T> & iv
+)
+{
+
+	str.readBegin("triple<T>");
+
+    str >> iv.x_;
+    str >> iv.y_;
+    str >> iv.z_;
+    
+
+    str.readEnd("triple<T>");
+
+    str.check(FUNCTION_NAME);
+
+	return str;
+}
+
+template<typename T>
+INLINE_FUNCTION void pFlow::readIstream
+(
+	iIstream & str,
+	triple<T> & iv
+)
+{
+
+	str.readBegin("triple<T>");
+	T val;
+
+    readIstream(str, val);
+    iv.x_ = val;
+    
+    readIstream(str, val);
+    iv.y_ = val;
+
+    readIstream(str, val);
+    iv.z_ = val;
+    
+    str.readEnd("triple<T>");
+
+    str.check(FUNCTION_NAME);
+
+}
+
+template<typename T>
+bool INLINE_FUNCTION_HD pFlow::equal
+(
+	const triple<T>& opr1,
+	const triple<T>& opr2
+)
+{
+	return equal( opr1.x(), opr2.x() ) && equal( opr1.y(), opr2.y() ) && equal( opr1.z(), opr2.z() );
+}
\ No newline at end of file
diff --git a/src/phasicFlow/types/triple/tripleMath.H b/src/phasicFlow/types/triple/tripleMath.H
new file mode 100644
index 00000000..00566ed3
--- /dev/null
+++ b/src/phasicFlow/types/triple/tripleMath.H
@@ -0,0 +1,93 @@
+/*------------------------------- 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.
+
+-----------------------------------------------------------------------------*/
+
+#define T3Func(fnName) 									\
+template<typename T> 									\
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::fnName(const pFlow::triple<T>& v) 			\
+{																\
+	return pFlow::triple<T>(fnName(v.x_), fnName(v.y_), fnName(v.z_));	\
+}
+
+#define T3Func2(fnName)  															\
+template<typename T>																\
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::fnName(const pFlow::triple<T>& arg1, const pFlow::triple<T>& arg2)	\
+{																					\
+	return pFlow::triple<T>(fnName(arg1.x_, arg2.x_), fnName(arg1.y_,arg2.y_), fnName(arg1.z_, arg2.z_)); \
+}
+
+//* * * * * * * * * * *  List of functinos * * * * * * * * //
+// abs, mod, exp, log, log10, pow, sqrt, cbrt
+// sin, cos, tan, asin, acos, atan, atan2
+// sinh, cosh, tanh, asinh, acosh, atanh 
+// min, max
+//* * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+T3Func(abs);
+T3Func2(mod);
+T3Func(exp);
+T3Func(log);
+T3Func(log10);
+T3Func2(pow);
+T3Func(sqrt);
+T3Func(cbrt);
+T3Func(sin);
+T3Func(cos);
+T3Func(tan);
+T3Func(asin);
+T3Func(acos);
+T3Func(atan);
+T3Func2(atan2);
+T3Func(sinh);
+T3Func(cosh);
+T3Func(tanh);
+T3Func(asinh);
+T3Func(acosh);
+T3Func(atanh);
+T3Func2(min);
+T3Func2(max);
+
+//// special forms of functions 
+
+// elements of t3 raised by e
+template<typename T>
+INLINE_FUNCTION_HD pFlow::triple<T> pFlow::pow(const pFlow::triple<T>& t3, T e)
+{
+	return pFlow::triple<T>( pow(t3.x_, e), pow(t3.y_,e), pow(t3.z_,e));
+}
+
+// return the min of 3 elements x, y, z
+template<typename T>
+INLINE_FUNCTION_HD T pFlow::min(const pFlow::triple<T>& t3)
+{
+	return min( min(t3.x_, t3.y_), t3.z_);
+}
+
+// return the max of 3 elements x, y, z
+template<typename T>
+INLINE_FUNCTION_HD T pFlow::max(const pFlow::triple<T>& t3)
+{
+	return max( max(t3.x_, t3.y_), t3.z_);
+}
+
+
+
+#undef T3Func
+#undef T3Func2
\ No newline at end of file
diff --git a/src/phasicFlow/types/types.C b/src/phasicFlow/types/types.C
new file mode 100644
index 00000000..6ebef81b
--- /dev/null
+++ b/src/phasicFlow/types/types.C
@@ -0,0 +1,39 @@
+/*------------------------------- 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 "types.H"
+
+namespace pFlow
+{
+		
+	const realx3 zero3(0.0);
+	const realx3 one3(1.0);
+	const uint32x3 zeroU3(0);
+	const uint32x3 oneU3(1);
+
+	const realx3x3   zero33(zero3);
+	const realx3x3   one33(one3);
+	const uint32x3x3 zeroU33(zeroU3);
+	const uint32x3x3 oneU33(oneU3);
+
+	const real4 zero4(zero);
+	
+
+} // pFlow
\ No newline at end of file
diff --git a/src/phasicFlow/types/types.H b/src/phasicFlow/types/types.H
new file mode 100644
index 00000000..e6d95268
--- /dev/null
+++ b/src/phasicFlow/types/types.H
@@ -0,0 +1,112 @@
+/*------------------------------- 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 __types_H__
+#define __types_H__
+
+
+#include "bTypes.H"
+
+#include "bTypesFunctions.H"
+
+#include "triple.H"
+
+#include "quadruple.H"
+
+#include "typeInfo.H"
+
+
+namespace pFlow
+{
+
+using int8x3 	= triple<int8>;
+using int16x3 	= triple<int16>;
+using int32x3 	= triple<int32>;
+using int64x3	= triple<int64>;
+using uint16x3 	= triple<uint16>;
+using uint32x3 	= triple<uint32>;
+using int32x3 	= triple<int32>;
+using int64x3 	= triple<int64>;
+using labelx3 	= triple<label>;
+using realx3 	= triple<real>;
+
+using uint16x3x3= triple<uint16x3>;
+using uint32x3x3= triple<uint32x3>;
+using int32x3x3 = triple<int32x3>;
+using labelx3x3 = triple<labelx3>;
+using realx3x3 	= triple<realx3>;
+
+using real4 	= quadruple<real>;
+
+
+template<>
+inline word basicTypeName<int8x3>(){	return "int8x3"; }
+
+template<>
+inline word basicTypeName<int16x3>(){	return "int16x3"; }
+
+template<>
+inline word basicTypeName<int32x3>(){	return "int32x3"; }
+
+template<>
+inline word basicTypeName<int64x3>(){	return "int64x3"; }
+
+template<>
+inline word basicTypeName<uint16x3>(){	return "uint16x3"; }
+
+template<>
+inline word basicTypeName<uint32x3>(){	return "uint32x3"; }
+
+template<>
+inline word basicTypeName<labelx3>(){	return "labelx3"; }
+
+template<>
+inline word basicTypeName<realx3>(){	return "realx3";	}
+
+template<>
+inline word basicTypeName<uint16x3x3>(){ return "uint16x3x3"; }
+
+template<>
+inline word basicTypeName<uint32x3x3>(){ return "uint32x3x3"; }
+
+template<>
+inline word basicTypeName<realx3x3>(){ return "realx3x3"; }
+
+
+template<>
+inline word basicTypeName<real4>(){	return "real4";	}
+
+
+extern const realx3 zero3;
+extern const realx3 one3;
+extern const uint32x3 zeroU3;
+extern const uint32x3 oneU3;
+
+extern const realx3x3 zero33;
+extern const realx3x3 one33;
+extern const uint32x3x3 zeroU33;
+extern const uint32x3x3 oneU33;
+
+
+
+} // pFlow
+
+
+#endif //__types_H__
\ No newline at end of file