new-message-observer-subscriber and anyList for holding a list of variables of any type

src/phasicFlow/CMakeLists.txt

@@ -20,7 +20,8 @@ streams/TStream/oTstream.cpp
 streams/Fstream/iFstream.cpp
 streams/Fstream/oFstream.cpp
 streams/Fstream/fileStream.cpp
-streams/dataIO/dataIO.cpp 
+streams/dataIO/dataIO.cpp
+streams/dataIO/dataIONoMPI.cpp 
 streams/streams.cpp 
 
 fileSystem/fileSystem.cpp
@@ -34,15 +35,22 @@ dictionary/twoPartEntry/twoPartEntry.cpp
 
 containers/Vector/Vectors.cpp
 containers/Field/Fields.cpp
+containers/List/anyList/anyList.cpp
 
 repository/IOobject/objectFile.cpp
 repository/IOobject/IOfileHeader.cpp
 repository/IOobject/IOobject.cpp
 repository/IOobject/IOPattern.cpp
 
+eventManagement/subscriber.cpp
+eventManagement/observer.cpp
+
 structuredData/line/line.cpp
 structuredData/infinitePlane/infinitePlane.cpp
 structuredData/plane/plane.cpp
+structuredData/domain/domain.cpp
+structuredData/pointStructure/internalPoints.cpp
+structuredData/boundaries/boundaryBase/boundaryBase.cpp
 
 )
 

src/phasicFlow/containers/Field/Field.cpp

@@ -216,10 +216,14 @@ bool pFlow::Field<VectorField, T, PropType>::write(
 	IOPattern::IOType iotype)const
 {
 
-	
 	os.writeWordKeyword(fieldKey_)<<endl;
-	VectorType::write(os, iotype);
+
+	if(!os.check(FUNCTION_NAME))return false;
+
+	if(!VectorType::write(os, iotype)) return false;
+
 	os.endEntry();
+	if(!os.check(FUNCTION_NAME))return false;
 
 	return true;
 }

src/phasicFlow/containers/Field/Field.hpp

@@ -197,6 +197,26 @@ public:
 		
 
 		bool write(iOstream& os, IOPattern::IOType iotype )const;
+
+		template<typename HostMask>
+		bool write(
+			iOstream& os, 
+			IOPattern::IOType iotype, 
+			const HostMask& mask)const
+		{
+			
+			os.writeWordKeyword(fieldKey_)<<endl;
+
+			if(!os.check(FUNCTION_NAME))return false;
+
+			if(!VectorType::write(os, iotype)) return false;
+
+			os.endEntry();
+			if(!os.check(FUNCTION_NAME))return false;
+
+			return true;
+		}
+		
 			
 
 };
@@ -205,7 +225,7 @@ public:
 template<template<class, class> class VectorField, class T, class PropType>
 inline iIstream& operator >> (iIstream & is, Field<VectorField, T, PropType> & ifld )
 {
-	if( !ifld.read(is, IOPattern::MasterProcessor) )
+	if( !ifld.read(is, IOPattern::MasterProcessorOnly) )
 	{
 		ioErrorInFile (is.name(), is.lineNumber());
 		fatalExit;

src/phasicFlow/containers/Field/Fields.hpp

@@ -34,7 +34,7 @@ using int8Field_D 		= Field<VectorSingle, int8>;
 
 using int8Field_H 		= Field<VectorSingle, int8, HostSpace>;
 
-/*using int32Field_D 		= Field<VectorSingle, int32>;
+using int32Field_D 		= Field<VectorSingle, int32>;
 
 using int32Field_H 		= Field<VectorSingle, int32, HostSpace>;
 
@@ -42,13 +42,17 @@ using int64Field_D 		= Field<VectorSingle, int64>;
 
 using int64Field_H 		= Field<VectorSingle, int64, HostSpace>;
 
+using uint8Field_D     = Field<VectorSingle, uint8>;
+
+using uint8Field_H     = Field<VectorSingle, uint8, HostSpace>;
+
 using uint32Field_D 	= Field<VectorSingle, uint32>;
 
 using uint32Field_H 	= Field<VectorSingle, uint32, HostSpace>;
 
-using labelField_D 		= Field<VectorSingle, label>;
+using uint64Field_D 		= Field<VectorSingle, uint64>;
 
-using labelField_H 		= Field<VectorSingle, label, HostSpace> ;
+using uint64Field_H 		= Field<VectorSingle, uint64, HostSpace> ;
 
 using realField_D 		= Field<VectorSingle, real>;
 
@@ -58,25 +62,9 @@ using realx3Field_D 	= Field<VectorSingle, realx3>;
 
 using realx3Field_H 	= Field<VectorSingle, realx3, HostSpace>;
 
-using uint16x3Field_D 	= Field<VectorSingle, uint16x3>;
-
-using uint16x3Field_H 	= Field<VectorSingle, uint16x3, HostSpace>;
-
-using uint32x3Field_D   = Field<VectorSingle, uint32x3>;
-
-using uint32x3Field_H   = Field<VectorSingle, uint32x3, HostSpace>;
-
-using int32x3Field_D 	= Field<VectorSingle, int32x3>;
-
-using int32x3Field_H 	= Field<VectorSingle, int32x3, HostSpace>;
-
-using int64x3Field_D 	= Field<VectorSingle, int64x3>;
-
-using int64x3Field_H 	= Field<VectorSingle, int64x3, HostSpace>;
-
 using realx3x3Field_D 	= Field<VectorSingle, realx3x3>;
 
-using realx3x3Field_H 	= Field<VectorSingle, realx3x3, HostSpace>;*/
+using realx3x3Field_H 	= Field<VectorSingle, realx3x3, HostSpace>;
 
 // - no typedef on device (since word does not compile on CUDA)
 using wordField_H 		= Field<VectorSingle, word, HostSpace>;

src/phasicFlow/containers/List/List/List.hpp

@@ -49,11 +49,11 @@ public:
 
 	using iterator 		= typename listType::iterator;
 	
-	using constIterator = typename listType::const_iterator;
+	using const_iterator = typename listType::const_iterator;
 	
 	using reference 	= typename listType::reference;
 	
-	using constReference= typename listType::const_reference;
+	using const_reference= typename listType::const_reference;
 	
 	using initList 		= typename std::initializer_list<T>;
 	
@@ -62,12 +62,7 @@ public:
 
 protected:
 
-	// position of ith element
+	
-	auto pos(size_t i);
-
-	// position of ith element 
-	const auto pos(size_t i)const;
-
 	static inline size_t getListStride(const size_t& len)
 	{
 		size_t stride = 1;
@@ -114,30 +109,19 @@ public:
 
 
 		// - copy construct
-		List(const List& src):
+		List(const List& src) = default;
-			listType(src)
+			
-		{}
 
 		// - move construct
-		List( List && mv)
+		List( List && mv) = default;
-		:
+		
-			listType(std::move(mv))
-		{}
-
 		// - copy assignment
-		ListType& operator=(const ListType& rhs)
+		ListType& operator=(const ListType& rhs) = default;
-		{
+		
-			listType::operator=(rhs);
-			return *this;
-		}
 
 		// - move assignment
-		ListType& operator=(ListType&& rhs)
+		ListType& operator=(ListType&& rhs) = default;
-		{
+		
-			listType::operator=(std::move(rhs));
-			return *this;
-		}	
-
 
 		uniquePtr<ListType> clone()const{
 			return makeUnique<ListType>(*this);
@@ -162,6 +146,12 @@ public:
 		// - size of container
 		size_t size()const;
 		
+		// position of ith element
+		auto pos(size_t i, bool noError = false);
+
+		// position of ith element 
+		const auto pos(size_t i, bool noError = false)const;
+
 		// - access to ith element
 		//   fatal exit if out of range
 		T& operator[](size_t i);
@@ -172,7 +162,7 @@ public:
 		
 		// - find the position of the first element with value val
 		//   cend() if not found
-		constIterator find(const T& val) const;
+		const_iterator find(const T& val) const;
 		
 		// - find the position of the first element with value val
 		//    end() if not found

src/phasicFlow/containers/List/List/ListI.hpp

@@ -21,14 +21,22 @@ Licence:
 template<typename T>
 auto pFlow::List<T>::pos
 (
-	size_t i
+	size_t i,
+	bool noError
 )
 {
 	if( i >= size() )
 	{
-		fatalErrorInFunction<<
+		if(noError)
-		"our of range access to list element. \n";
+		{
-		fatalExit;
+			return std::end(*this);
+		}
+		else
+		{
+			fatalErrorInFunction<<
+			"out of range access to list element. \n";
+			fatalExit;
+		}
 	}
 	auto iter = listType::begin();
 	std::advance(iter, i);
@@ -38,14 +46,22 @@ auto pFlow::List<T>::pos
 template<typename T>
 const auto pFlow::List<T>::pos
 (
-	size_t i
+	size_t i,
+	bool noError
 )const 
 {
 	if( i >= size() )
 	{
-		fatalErrorInFunction<<
+		if(noError)
-		"our of range access to list element. \n";
+		{
-		fatalExit;
+			return std::end(*this);
+		}
+		else
+		{
+			fatalErrorInFunction<<
+			"out of range access to list element. \n";
+			fatalExit;
+		}
 	}
 	auto iter = listType::cbegin();
 	std::advance(iter, i);
@@ -88,7 +104,7 @@ inline const T& pFlow::List<T>::operator[]
 }
 
 template<typename T>
-inline typename  pFlow::List<T>::constIterator pFlow::List<T>::find
+inline typename  pFlow::List<T>::const_iterator pFlow::List<T>::find
 (
 	const T& val
 ) const

src/phasicFlow/containers/List/anyList/anyList.cpp

@@ -0,0 +1,43 @@
+/*------------------------------- 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 "anyList.hpp"
+
+
+size_t pFlow::anyList::size()const
+{
+	return anyList_.size();
+}
+
+bool pFlow::anyList::empty()const
+{
+	return size()==0;
+}
+
+bool pFlow::anyList::contains(const word& name)const
+{
+	return names_.search(name);
+}
+
+const pFlow::wordList& pFlow::anyList::names()const
+{
+	return names_;
+}

src/phasicFlow/containers/List/anyList/anyList.hpp

@@ -0,0 +1,243 @@
+/*------------------------------- 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 __anyList_hpp__
+#define __anyList_hpp__
+
+#include <any>
+
+#include "types.hpp"
+#include "List.hpp"
+
+
+namespace pFlow
+{
+
+
+class anyList
+{
+
+public:
+
+	using anyListType 	= List<std::any>;
+
+	using reference 	= typename anyListType::reference;
+
+	using iterator 		= typename anyListType::iterator;
+
+	using const_iterator= typename anyListType::const_iterator;
+
+
+protected:
+
+	/// Contains a list of variables with any type 
+	anyListType 	anyList_;
+
+	/// List of variable names in anyList_ 
+	wordList 		names_;	
+
+public:
+
+	// - Type info
+	TypeInfoNV("anyList");
+
+	//// - Constructors
+
+		// All five constructors are created by compiler 
+	
+
+	//// -  Methods
+		
+		/// Size of container
+		size_t size()const;
+
+		/// Is container empty
+		bool empty()const;
+
+		/// Does container contain this variable name?
+		bool contains(const word& name)const;
+
+		/// List of varibales names
+		const wordList& names()const;
+
+		/// Create variable using constructor in-place 
+		template<typename T, typename... Args>
+		reference emplaceBack(const word& name, Args&&... args)
+		{
+			if( contains(name))
+			{
+				fatalErrorInFunction<<
+				"variable name "<< name << " already exists in the anyList."<<endl<<
+				"list of variables is "<<names_<<endl;
+				fatalExit;
+			}
+			names_.push_back(name);
+			return anyList_.emplace_back( 
+				std::in_place_type<T>, 
+				std::forward<Args>(args)...);
+		}
+
+		/// Create variable using copy constructor
+		template<typename T>
+		reference emplaceBack(const word& name, const T & other)
+		{
+			if( contains(name) )
+			{
+				fatalErrorInFunction<<
+				"variable name "<< name << " already exists in the anyList."<<endl<<
+				"list of variables is "<<names_<<endl;
+				fatalExit;
+			}
+			names_.push_back(name);
+			return anyList_.emplace_back(std::in_place_type<T>, other);
+		}
+
+		/// Create variable using move constructor 
+		template<typename T>
+		reference emplaceBack(const word& name, T&& other)
+		{
+			if( contains(name) )
+			{
+				fatalErrorInFunction<<
+				"variable name "<< name << " already exists in the anyList."<<endl<<
+				"list of variables is "<<names_<<endl;
+				fatalExit;
+			}
+			names_.push_back(name);
+			return anyList_.emplace_back(std::in_place_type<T>, std::move(other));
+		}
+
+		/// Get the reference to variable by index
+		template<typename T>
+		T& getObject(size_t i)
+		{
+			std::any *a = &(*anyList_.pos(i));
+			if(!a->has_value())
+			{
+				fatalErrorInFunction<<
+				"any does not have a value for dereferencing. "<<
+				"index in anyList is "<<i<<"."<<endl;
+				fatalExit;
+			}
+			return *std::any_cast<T>(a);
+		}
+
+		/// Get the reference to variable by name 
+		template<typename T>
+		T& getObject(const word& name)
+		{
+			int32 i = names_.findi(name);
+			if(i == -1 )
+			{
+				fatalErrorInFunction<<
+				"variable name "<< name << " does not exist in the anyList."<<endl<<
+				"list of variables is "<<names_<<endl;
+				fatalExit;
+			}
+			std::any *a = &(*anyList_.pos(i));
+			if(!a->has_value())
+			{
+				fatalErrorInFunction<<
+				"any does not have a value for dereferencing. "<<
+				"index in anyList is "<<i<<"."<<endl;
+				fatalExit;
+			}
+			return *std::any_cast<T>(a);
+		}
+
+		/// Get the const reference to variable by name 
+		template<typename T>
+		const T& getObject(const word& name)const
+		{
+			int32 i = names_.findi(name);
+			if(i == -1 )
+			{
+				fatalErrorInFunction<<
+				"variable name "<< name << " does not exist in the anyList."<<endl<<
+				"list of variables is "<<names_<<endl;
+				fatalExit;
+			}
+			const std::any *a = &(*anyList_.pos(i));
+			if(!a->has_value())
+			{
+				fatalErrorInFunction<<
+				"any does not have a value for dereferencing. "<<
+				"index in anyList is "<<i<<"."<<endl;
+				fatalExit;
+			}
+			return *std::any_cast<const T>(a);
+		}
+
+		/// Get the const reference to object by index
+		template<typename T>
+		const T& getObject(size_t i)const
+		{
+			const std::any *a = &(*anyList_.pos(i));
+			if(!a->has_value())
+			{
+				fatalErrorInFunction<<
+				"any does not have a value for dereferencing. "<<
+				"index in anyList is "<<i<<"."<<endl;
+				fatalExit;
+			}
+			return *std::any_cast<const T>(a);
+		}
+
+		/// Get the pointer to variable by index
+		template<typename T>
+		T* getObjectPtr(size_t i)
+		{
+			if(i>= size())return nullptr;
+
+			std::any *a = &(*anyList_.pos(i));
+			if( a->has_value())
+			{
+				return std::any_cast<T>(a);	
+			}
+			else
+			{
+				return nullptr;
+			}	
+		}
+
+		/// Get the const pointer to variable by index 
+		template<typename T>
+		const T* getObjectPtr(size_t i)const
+		{
+			if(i>=size())return nullptr;
+			const std::any *a = &(*anyList_.pos(i));
+			if( a->has_value())
+			{
+				return std::any_cast<const T>(a);	
+			}
+			else
+			{
+				return nullptr;
+			}	
+		}
+			
+};
+
+
+
+} // pFlow
+
+
+#endif //__anyList_hpp__

src/phasicFlow/containers/Vector/Vector.hpp

@@ -364,7 +364,7 @@ public:
 template<typename T, typename Allocator>
 inline iIstream& operator >> (iIstream & is, Vector<T, Allocator> & ivec )
 {
-	if( !ivec.readVector(is, IOPattern::MasterProcessor) )
+	if( !ivec.readVector(is, IOPattern::MasterProcessorOnly) )
 	{
 		ioErrorInFile (is.name(), is.lineNumber());
 		fatalExit;

src/phasicFlow/containers/Vector/stdVectorHelper.hpp

@@ -57,9 +57,15 @@ bool writeSpan(
     IOPattern::IOType iotype)
 {
    
-   dataIO io(iotype);
+   auto ioPtr = dataIO::create(iotype, IOPattern::exeMode());
 
-   if(!io.writeData(os, sp))
+   if(!ioPtr)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+   if(!ioPtr().writeData(os, sp))
    {
    		fatalErrorInFunction;
    		return false;
@@ -92,9 +98,15 @@ bool readStdVector
 	IOPattern::IOType iotype 
 )
 {
-	dataIO io(iotype);
+	auto ioPtr = dataIO::create(iotype, IOPattern::exeMode());
 
-	if(!io.readData(is, vec))
+	if(!ioPtr)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	if(!ioPtr().readData(is, vec))
 	{
 		fatalErrorInFunction;
 		return false;
@@ -118,7 +130,7 @@ iOstream& operator<<( iOstream& os, const std::vector<T,Allocator>& vec)
 template<typename T, typename Allocator>
 iIstream& operator>>(iIstream& is, std::vector<T,Allocator>& vec)
 {
-	if( !readStdVector(is,vec, IOPattern::MasterProcessor))
+	if( !readStdVector(is,vec, IOPattern::MasterProcessorOnly))
 	{
 		fatalErrorInFunction;
 		fatalExit;

src/phasicFlow/containers/VectorHD/VectorSingle.hpp

@@ -599,12 +599,38 @@ public:
 			
 		}
 
+		template<typename HostMask>
+		FUNCTION_H
+		bool write(iOstream& os, IOPattern::IOType iotype, const HostMask& mask)const
+		{
+			auto hVec = hostVector();
+			
+			auto numActive = mask.numActive();
+			std::vector<T> finalField;
+			finalField.clear();
+			finalField.reserve(numActive);
+			
+			uint32 start 	= mask.activeRange().start();
+			uint32 end 		= mask.activeRange().end();
+			
+			for(uint32 i=start; i<end; i++ )
+			{
+				if( mask() )
+					finalField.push_back(hVec[i]);
+			}
+
+			auto sp = span<T>( finalField.data(), finalField.size());
+
+			return writeSpan(os, sp, iotype);
+
+		}
+
 }; // class VectorSingle
 
 template<typename T, typename MemorySpace>
 inline iIstream& operator >> (iIstream & is, VectorSingle<T, MemorySpace> & ivec )
 {
-	if( !ivec.read(is, IOPattern::MasterProcessor ) )
+	if( !ivec.read(is, IOPattern::MasterProcessorOnly ) )
 	{
 		ioErrorInFile (is.name(), is.lineNumber());
 		fatalExit;

src/phasicFlow/containers/span/span.hpp

@@ -77,11 +77,11 @@ public:
 
     /// move
     INLINE_FUNCTION_HD
-    span(span&&) = delete;
+    span(span&&) = default;
 
     /// assignment
     INLINE_FUNCTION_HD
-    span& operator=(span&) = delete;
+    span& operator=(span&) = default;
 
     
     INLINE_FUNCTION_HD

src/phasicFlow/eventManagement/message.hpp

@@ -0,0 +1,159 @@
+/*------------------------------- 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 __message_hpp__
+#define __message_hpp__
+
+#include <bitset>
+
+#include "types.hpp"
+#include "iOstream.hpp"
+
+
+namespace pFlow
+{
+
+class message
+{
+public:
+	enum EVENT : size_t
+	{
+		DEFAULT 		= 0,
+		CAP_CHANGED 	= 1,
+		SIZE_CHANGED 	= 2,
+		ITEM_DELETE	 	= 3,
+		ITEM_INSERT  	= 4,
+		RANGE_CHANGED   = 5,
+		ITEM_REARRANGE 	= 6	
+	};
+
+protected:
+
+	std::bitset<16> events_{0x0000};
+	
+
+public:
+
+	message()=default;
+
+	message(EVENT evnt)
+	{
+		add(evnt);
+	}
+
+	message(size_t i )
+	{
+		if(i<16)
+		{
+			events_.set(i);
+		}
+	}
+
+	message(const message&) = default;
+
+	message(message&&) = default;
+
+	message& operator=(const message&) = default;
+
+	message& operator=(message&&) = default;
+
+	~message()=default;
+
+	inline message& add( EVENT evnt)
+	{
+		events_.set(static_cast<size_t>(evnt));
+		return *this;
+	}
+
+	inline message& remove(EVENT evnt)
+	{
+		events_.set(static_cast<size_t>(evnt), false);
+		return *this;
+	}
+
+	inline bool equivalentTo( EVENT evnt )const
+	{
+		return equivalentTo(static_cast<size_t>(evnt));
+	}
+
+	inline bool equivalentTo(size_t i)const
+	{
+		return events_.test(i);
+	}
+
+	inline bool isNull()const
+	{
+		return events_.none();
+	}
+
+	inline auto size()const
+	{
+		return events_.size();
+	}
+
+	inline
+	message& operator+(EVENT evnt)
+	{
+		return add(evnt);
+	}
+
+	inline
+	message& operator-(EVENT evnt)
+	{
+		return remove(evnt);
+	}
+
+	static 
+	message Default()
+	{
+		message msg;
+		return msg+DEFAULT;
+	}
+
+	static 
+	message Empty()	
+	{
+		message msg;
+		return msg;
+	}
+	
+};
+
+inline
+iOstream& operator<<(iOstream& os, const message& msg)
+{
+	for(size_t i=msg.size(); i>0; i--)
+	{
+		os<< msg.equivalentTo(i-1)? '1':0;
+	}
+	os<<endl;
+
+	return os;
+}
+
+inline
+message operator + (message::EVENT evnt1, message::EVENT evnt2)
+{
+	message msg = message::Empty().add(evnt1).add(evnt2);
+	return msg;
+}
+
+}
+
+#endif // __eventMessage_hpp__

src/phasicFlow/eventManagement/observer.cpp

@@ -0,0 +1,59 @@
+/*------------------------------- 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 "observer.hpp"
+#include "subscriber.hpp" 
+
+pFlow::observer::observer():
+	subscriber_(nullptr)
+{}
+
+pFlow::observer::observer
+(
+	const subscriber* subscrbr,
+	message msg
+)
+:
+	subscriber_(subscrbr),
+	message_(msg)
+{
+	if(subscriber_)
+	{
+		if(!subscriber_->subscribe(msg, this))
+		{
+			fatalErrorInFunction<<
+			"error in subscribing an observer"<<endl;
+			fatalExit;
+		}
+	}
+}
+
+pFlow::observer::~observer()
+{
+	if( subscriber_)
+		subscriber_->unsubscribe(this);
+	invalidateSubscriber();
+}
+
+bool pFlow::observer::addToSubscriber(const subscriber& subscrbr)
+{
+	subscriber_ = &subscrbr;
+	return subscriber_->subscribe(message_, this);
+}

src/phasicFlow/eventManagement/observer.hpp

@@ -0,0 +1,76 @@
+/*------------------------------- 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 __observer_hpp__
+#define __observer_hpp__
+
+#include "message.hpp"
+
+namespace pFlow
+{
+
+class subscriber;
+class anyList;
+
+class observer
+{
+protected:
+	
+
+	/// pointer to subscriber 
+	const subscriber* 	subscriber_ = nullptr;
+	
+	/// list of events in the message 
+	const message 		message_;
+
+
+public:
+
+	observer();
+
+	observer(
+		const subscriber* subscrbr, 
+		message msg);
+
+	virtual 
+	~observer();
+
+	inline 
+	bool subscribed()const 
+	{
+		return subscriber_!=nullptr;
+	}
+
+	bool addToSubscriber(const subscriber& subscriber);
+	
+	inline void invalidateSubscriber()
+	{
+		subscriber_ = nullptr;
+	}
+
+
+	virtual bool update(const message& msg, const anyList& varList)=0;
+};
+
+} // pFlow
+
+
+#endif // __observer_hpp__

src/phasicFlow/eventManagement/subscriber.cpp

@@ -0,0 +1,147 @@
+/*------------------------------- 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 "subscriber.hpp"
+#include "anyList.hpp"
+#include "observer.hpp"
+#include "message.hpp"
+
+
+pFlow::subscriber::~subscriber()
+{
+	for(size_t i=0; i<observerList_.size(); i++)
+	{
+		for( auto obsvr: observerList_[i] )
+		{
+			if(obsvr)
+			{
+				obsvr->invalidateSubscriber();
+			}
+		}
+	}
+}
+
+bool pFlow::subscriber::subscribe
+(
+	message 	msg,
+	observer* 	obsevr
+)const
+{
+	if( msg.size()>16 )
+	{
+		fatalErrorInFunction<<
+		"message size is greater than 16!"<<endl;
+		return false;
+	}
+	// do not subscribe nullptr 
+	if(!obsevr)return false;
+
+	for(size_t i=0; i<msg.size(); i++)
+	{
+		if(msg.equivalentTo(i))
+		{
+			if(!observerList_[i].search(obsevr))
+			{
+				observerList_[i].push_back(obsevr);
+			}
+		}
+	}
+	return true;
+}
+
+bool pFlow::subscriber::unsubscribe
+(
+	observer* obsrvr
+)const
+{
+	if(obsrvr)
+	{
+		for(size_t i=0; i<observerList_.size(); i++)
+		{
+			if( auto iter = observerList_[i].find(obsrvr); iter!= observerList_[i].end())
+			{
+				(*iter)->invalidateSubscriber();
+				observerList_[i].erase(iter);
+			}
+		}
+	}
+	return true;
+}
+
+bool pFlow::subscriber::notify
+(
+	const message msg, 
+	const anyList& varList
+)
+{
+	if( msg.size()>16 )
+	{
+		fatalErrorInFunction<<
+		"message size is greater than 16!"<<endl;
+		return false;
+	}
+
+	for(size_t i=0; i<msg.size(); i++)
+	{
+		if(msg.equivalentTo(i))
+		{
+			for( auto obsvr: observerList_[i] )
+			{
+				obsvr->update(message(i), varList);
+			}
+		}
+	}
+
+	return true;
+}
+
+/*bool pFlow::subscriber::notify
+(
+	const eventMessage &msg
+)
+{
+	for ( auto& observer:observerList_ )
+	{
+		if(observer)
+			if( !observer->update(msg) ) return false;
+	}
+
+	return true;
+}
+
+bool pFlow::eventSubscriber::notify
+(
+	const eventMessage& msg,
+	const List<eventObserver*>& exclutionList
+)
+{
+	Set<eventObserver*> sortedExcList(exclutionList.begin(),exclutionList.end());
+
+	for(auto& observer:observerList_)
+	{
+		if( observer && sortedExcList.count(observer) == 0 )
+		{
+			if(!observer->update(msg)) return false;
+		}
+	}
+
+	return true;
+}*/

src/phasicFlow/eventManagement/subscriber.hpp

@@ -0,0 +1,68 @@
+/*------------------------------- 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 __subscriber_hpp__
+#define __subscriber_hpp__
+
+
+
+#include <array>
+
+#include "List.hpp"
+
+namespace pFlow
+{
+
+class observer;
+class message;
+class anyList;
+
+class subscriber
+{
+protected:
+
+	// - list of subsribed objectd that recieve updage messages 
+	mutable std::array<List<observer*>,16> observerList_;
+
+public:
+
+	subscriber()
+	{}
+
+	virtual ~subscriber();
+
+	virtual bool subscribe(message msg, observer* obsevr)const;
+
+	virtual bool unsubscribe(observer* obsevr)const;
+
+	//bool notify(const eventMessage& msg);
+
+	//bool notify(const eventMessage& msg, const List<eventObserver*>& exclutionList );
+
+	
+	bool notify(const message msg, const anyList& varList);
+	
+
+};
+
+} // pFlow
+
+
+
+#endif // __eventSubscriber_hpp__

src/phasicFlow/repository/IOobject/IOPattern.cpp

@@ -29,3 +29,10 @@ pFlow::IOPattern::IOPattern( IOType iotype)
 	isParallel_(processors::isParallel())
 {}
 
+pFlow::word pFlow::IOPattern::exeMode()
+{
+    if(processors::isParallel())
+    	return word("MPI");
+    else
+    	return word("NoMPI");
+}

src/phasicFlow/repository/IOobject/IOPattern.hpp

@@ -22,7 +22,8 @@ Licence:
 #define __IOPattern_hpp__
 
 
-#include "processors.hpp"
+
+#include "types.hpp"
 
 namespace pFlow
 {
@@ -52,7 +53,9 @@ public:
 	{
 		MasterProcessorOnly			= 0,
 		AllProcessorsSimilar 		= 1,
-		MasterProcessorDistribute	= 4
+		MasterProcessorDistribute	= 4,
+		AllProcessorsDifferent  = 8 // this is used for << and >> operators for 
+									// standard input and output streams 
 	};
 
 protected:
@@ -95,6 +98,12 @@ public:
 		return ioType_ == MasterProcessorDistribute;
 	}
 
+	inline
+	bool isAllProcessorsDifferent()const
+	{
+		return ioType_ == AllProcessorsDifferent;
+	}
+
 	inline
 	bool isMaster()const
 	{
@@ -110,7 +119,7 @@ public:
 	inline
 	bool thisProcReadData()const
 	{
-		if(isMasterProcessor() && !isMaster())return false;
+		if(isMasterProcessorOnly() && !isMaster())return false;
 		if(isMasterProcessorDistribute() && !isMaster()) return false;
 		return true;
 	}
@@ -118,6 +127,7 @@ public:
 	inline
 	bool thisProcWriteData()const
 	{
+		if(isAllProcessorsDifferent()) return true;
 		return isMaster();
 	}
 
@@ -145,6 +155,10 @@ public:
 		return globalRank_;
 	}
 
+	static 
+	word exeMode();
+
+
 };
 
 }

src/phasicFlow/repository/IOobject/objectFile.hpp

@@ -66,7 +66,7 @@ protected:
 	/// Number of bytes used for writing/reading real variable (used for binray)
 	int 		numBytesForReal_ = numBytesForReal__;
 	
-	IOPattern 	ioPattern_ = {IOPattern::MasterProcessor};
+	IOPattern 	ioPattern_ = {IOPattern::MasterProcessorOnly};
 
 	/// Does the objectFile write the header or not
 	bool readWriteHeader_ = true;
@@ -88,7 +88,7 @@ public:
 		const fileSystem& 	localPath,
 		const readFlag&   	rf = READ_NEVER,
 		const writeFlag&  	wf = WRITE_NEVER,
-		IOPattern::IOType	ioType = IOPattern::MasterProcessor,
+		IOPattern::IOType	ioType = IOPattern::MasterProcessorOnly,
 		bool  rwHdr = true
 	);
 

src/phasicFlow/smartPointers/uniquePtr.hpp

@@ -54,24 +54,20 @@ public:
 
 
 	template<typename... Args>
-    inline static uniquePtr<T> makeUnique(Args&&... args)
+	inline static uniquePtr<T> makeUnique(Args&&... args)
-    {
+	{
-    	return uniquePtr<T>(new T(std::forward<Args>(args)...));
+		return uniquePtr<T>(new T(std::forward<Args>(args)...));
-    }
+	}
 
-    void clear()
+	void clear()
-    {
+	{
-    	if( ! (*this) )
+		this->reset(nullptr);
-    	{
+	}
-    		auto p = this->release();
-    		delete p;
-    	}
-    }
 
-    // ref to the object of type T
+	// ref to the object of type T
 	T& operator ()()
 	{
-		if(!this->get())
+		if(!*this)
 		{
 			fatalErrorInFunction <<
 			"uniquePtr is empty, and you are trying to get its reference. \n" <<
@@ -84,14 +80,19 @@ public:
 	// const ref to the object of type T
 	const T& operator() () const
 	{
-		if(!this->get())
+		if(!*this)
 		{
 			fatalErrorInFunction <<
 			"uniquePtr is empty, and you are trying to get its reference. \n" <<
 			"Type name is "<< typeid(T).name()<<"\n";
 			fatalExit;
 		}
-		return const_cast<T&>(*this->get());
+		return static_cast<const T&>(*this->get());
+	}
+
+	explicit operator bool() const
+	{
+		return this->get()!= nullptr;
 	}
 
 };

src/phasicFlow/streams/dataIO/dataIO.cpp

@@ -1,717 +1,34 @@
 
-#include "phasicFlowConfig.H"
-#ifdef pFlow_Build_MPI
-	#include <mpi.h>
-#endif
-
-#include <cstdio>
-#include <numeric>
-
 #include "dataIO.hpp"
 
 
-
+pFlow::uniquePtr<pFlow::dataIO> pFlow::dataIO::create
-static const size_t numcharFlag = 8;
-
-static inline unsigned char binaryFlag__[numcharFlag] = {255, 255, 255, 255, 255, 255, 255 ,0};
-
-inline 
-pFlow::uint64 chunkSizeOffeset(pFlow::uint64 procNo)
-{
-	// 1 is added for number of processors 
-	return (1+procNo)*sizeof(pFlow::uint64);
-}
-
-bool writeBinaryBlockFlagSTD(std::FILE* fh)
-{
-
-	std::size_t nw = std::fwrite( binaryFlag__ , sizeof(unsigned char), numcharFlag, fh );
-
-	if(nw < numcharFlag) 
-		return false;
-	else
-		return true;
-}
-
-pFlow::uint64 findBindaryBlockFlagSTD(std::FILE* fh)
-{
-	// get the current postion 
-	long fpos;
-	if( fpos = std::ftell( fh) ; fpos == -1L )
-	{
-		fatalErrorInFunction;
-		return pFlow::dataIO::ErrorReturn;
-	}
-
-	pFlow::uint64 filePos = static_cast<pFlow::uint64>(fpos);
-
-	// start reading char by char 
-	unsigned char ch;
-	int currPos = 0;
-    while ( std::fread(&ch, sizeof(ch), 1, fh) == 1 ) 
-    {
-   		if(std::ferror(fh)) return pFlow::dataIO::ErrorReturn;	
-   		if(std::feof(fh))return pFlow::dataIO::ErrorReturn;
-
-   		filePos++;
-
-   		if(ch == binaryFlag__[currPos] ) 
-   		{
-   			currPos++;
-
-   			if(currPos == numcharFlag) return filePos;
-   		}
-   		else
-   		{
-   			currPos = 0;
-   		}   
-    }
-
-    return pFlow::dataIO::ErrorReturn;
-}
-
-
-#ifdef pFlow_Build_MPI
-pFlow::uint64 findBindaryBlockFlagMPI(MPI_File fh, pFlow::uint64 startPosSearch)
-{
-	
-
-	pFlow::uint64 filePos = static_cast<pFlow::uint64>(startPosSearch);
-
-	// start reading char by char 
-	unsigned char ch;
-	int currPos = 0;
-
-    while( MPI_File_read(
-    	fh, 
-    	&ch, 
-    	1, 
-    	MPI_UNSIGNED_CHAR, 
-    	MPI_STATUS_IGNORE) == MPI_SUCCESS ) 
-    {
-   		
-   		filePos++;
-   		
-   		if(ch == binaryFlag__[currPos] ) 
-   		{
-   			currPos++;
-
-   			if(currPos == numcharFlag) return filePos;
-   		}
-   		else
-   		{
-   			currPos = 0;
-   		}   
-    }
-
-    return pFlow::dataIO::ErrorReturn;
-}
-
-#endif
-
-
-bool pFlow::dataIO::writeDataToFileEndSTD
 (
-	const word& wordPath, 
+	IOPattern::IOType iotype, 
-	const span<unsigned char>& data
+	word exeMode
 )
 {
+	auto ioMode = angleBracketsNames("dataIO", exeMode);
 
-	// openfile 
+	if( wordvCtorSelector_.search(ioMode) )
-	auto fh = std::fopen(wordPath.c_str(), "ab");
-
-	if(!fh)
 	{
-		fatalErrorInFunction<<
+		auto objPtr = wordvCtorSelector_[ioMode] 
-		"Error in Opening file "<< wordPath <<endl;
+			(
-		std::fclose(fh);
+				iotype,
-		return false;
+				exeMode
-	}
+			);
-
+		return objPtr;
-	// write the start of binary block flag
-	if(std::fseek(fh, 0 , SEEK_END)!=0)
-	{
-		fatalErrorInFunction<<
-		"error at reaching end of file "<<wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-	
-	
-	if(!writeBinaryBlockFlagSTD(fh) )
-	{
-		fatalErrorInFunction<<
-		"Error in writing to file "<< wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-		
-	uint64 numChunks = 1;	
-
-	// write number of data chunks
-	auto wc = std::fwrite(&numChunks, sizeof(numChunks), 1, fh);
-
-	if(wc < 1)
-	{
-		fatalErrorInFunction<<
-		"Error in writing numChunks to file "<< wordPath <<endl;
-		std::fclose(fh);
-		return false;
-	}
-
-
-	// write size of each data chunk
-	uint64 sizeOfData = data.size();
-
-	wc = std::fwrite(&sizeOfData, sizeof(sizeOfData), 1, fh);
-	if(wc <1)
-	{
-		fatalErrorInFunction<<
-		"Error in writing size of data chunk to file "<< wordPath <<endl;
-		std::fclose(fh);
-		return false;	
-	}
-
-	if(sizeOfData>0 )
-	{
-		// write data chunk to file
-		wc = std::fwrite(data.data(), sizeof(unsigned char), sizeOfData, fh);
-		if(wc < sizeOfData )
-		{
-			fatalErrorInFunction<<
-			"Error in writing size of data to file "<< wordPath <<endl;
-			std::fclose(fh);
-			return false;	
-		}	
-	}
-	
-		
-	lastPosRead_ = static_cast<uint64>(std::ftell(fh));
-		
-	// close the file 
-	std::fclose(fh);
-
-	return true;
-}
-
-bool pFlow::dataIO::writeDataToFileEndMPI
-(
-	const word& wordPath, 
-	const span<unsigned char>& data
-)
-{	
-
-#ifdef pFlow_Build_MPI
-	
-	// collect information from other processes
-	uint64 numProc 	= processors::globalSize();
-	uint64 thisSize = data.size();
-	uint64 offset;
-
-	CheckMPI
-	(
-		MPI_Scan(&thisSize, &offset, 1, MPI_UINT64_T,  MPI_SUM, MPI_COMM_WORLD),
-		true
-	);
-
-
-	MPI_File fh;
-	
-	if( MPI_File_open(
-			MPI_COMM_WORLD,
-			wordPath.c_str(),
-			MPI_MODE_WRONLY+MPI_MODE_APPEND,
-			MPI_INFO_NULL,
-			&fh) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction<<
-		"Cannot open file "<< wordPath<<endl;
-		return false;
-	}	
-
-	uint64 startPos;
-
-	if( processors::isMaster())
-	{
-		// set position to the end of file 
-		MPI_File_seek(fh, 0 , MPI_SEEK_END);		
-
-		if( MPI_File_write(
-			fh, 
-			binaryFlag__,
-			numcharFlag,
-			MPI_UNSIGNED_CHAR,
-			MPI_STATUS_IGNORE
-			) != MPI_SUCCESS )
-		{
-			fatalErrorInFunction<<
-			"Cannot write binary block flag into "<< wordPath<<endl;
-			return false;
-		}
-		
-		MPI_Offset posOfBlock;
-		if( MPI_File_get_position(fh, &posOfBlock) != MPI_SUCCESS )
-		{
-			fatalErrorInFunction<<
-			"Cannot get the end pos of file "<< wordPath<<endl;
-			return false;
-		}
-
-		
-		startPos = static_cast<uint64>(posOfBlock);
-	}
-
-		
-	if( MPI_Bcast(
-		&startPos, 
-		1, 
-		MPI_UINT64_T, 
-		processors::masterNo(), 
-		MPI_COMM_WORLD) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction;
-		return false;
-	}
- 
-
-	if( processors::isMaster())
-	{
-		uint64 numProc = processors::globalSize();
-		if(MPI_File_write(
-			fh,
-			&numProc,
-			1,
-			MPI_UINT64_T,
-			MPI_STATUS_IGNORE) != MPI_SUCCESS)
-		{
-			fatalErrorInFunction<<
-			"Cannot write number of chunks into "<<wordPath<<endl;
-			return false;
-		}
-	}
-
-	uint64 sizeOffset = startPos  + chunkSizeOffeset(processors::globalRank());
-
-	
-	if(MPI_File_write_at_all(
-		fh,
-		sizeOffset, 
-		&thisSize,  
-		1,
-		MPI_UINT64_T,
-		MPI_STATUS_IGNORE)!= MPI_SUCCESS)
-	{
-		fatalErrorInFunction<<
-		"Cannot write size of chunk into "<<wordPath<<endl;
-		return false;
-	}
-
-	offset -= thisSize;
-
-	uint64 chunkOffset = startPos + chunkSizeOffeset(processors::globalSize()) + offset;
-
-	
-	if(MPI_File_write_at_all(
-		fh,
-		chunkOffset,
-		data.data(),
-		thisSize,
-		MPI_UNSIGNED_CHAR,
-		MPI_STATUS_IGNORE) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction<<
-		"Cannot write data into "<<wordPath<<endl;
-		return false;;
-	}
-
-	
-
-	MPI_File_close(&fh);
-	MPI_Barrier(MPI_COMM_WORLD);
-
-	uint64 lastPos = chunkOffset + thisSize;
-	return maxReduction(lastPos, lastPosWrite_);
-
-	/*if( MPI_Allreduce(
-			&lastPos,
-			&lastPosWrite_,
-			1,
-			MPI_UINT64_T,
-			MPI_MAX,
-			MPI_COMM_WORLD) != MPI_SUCCESS )
-	{
-		fatalErrorInFunction<<
-		"Error in max_reduction for write operation"<<endl;
-		return false;
-	}
-
-	return true;*/
-
-#else
-	// if MPI is not active, we use std routins
-	return writeDataToFileEndSTD(wordPath, data);
-#endif
-}
-
-
-bool pFlow::dataIO::readDataSTD
-(
-	const word& wordPath,
-	const std::vector<uint64> chunkSizes,
-	span<unsigned char>& data,
-	uint64 binaryBlockStart
-)
-{
-	
-	
-	// sum of all chuncks
-	uint64 toRecv = std::accumulate(
-		chunkSizes.begin(), 
-		chunkSizes.end(), 
-		static_cast<uint64>(0));
-
-	if( data.size() != toRecv )
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-
-	
-	auto fh = std::fopen(wordPath.c_str(), "rb");
-
-	if(!fh)
-	{
-		fatalErrorInFunction<<
-		"Error in Opening file "<< wordPath<<endl;
-		return false;
-	}
-
-	// start of data chunks 
-	uint64 offset = binaryBlockStart + chunkSizeOffeset(chunkSizes.size());
-	
-
-	if(auto res = std::fseek(fh, offset, SEEK_SET); res!= 0 )
-	{
-		fatalErrorInFunction<<
-		"Error in file seek "<< wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-
-
-	if(auto res = std::fread(
-			data.data(),
-			sizeof(unsigned char),
-			data.size(),
-			fh);
-		res!= data.size() )
-	{
-		fatalErrorInFunction<<
-		"Error in reading file "<< wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-
-	lastPosRead_ = static_cast<uint64>(std::ftell(fh));
-		
-	std::fclose(fh);
-	return true;
-}
-
-
-bool pFlow::dataIO::readDataMPI
-(
-	const word& wordPath, 
-	const std::vector<uint64> chunkSizes,
-	span<unsigned char>& data,
-	uint64 binaryBlockStart
-)
-{
-
-#ifdef pFlow_Build_MPI
-
-	
-	MPI_File fh;
-
-	if( MPI_File_open(
-		MPI_COMM_WORLD,
-		wordPath.c_str(),
-		MPI_MODE_RDONLY,
-		MPI_INFO_NULL,
-		&fh))
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-		
-	
-	auto procNo = processors::globalRank();
-
-	uint64 toRecv = chunkSizes[procNo];
-
-	// start of data chunks 
-	uint64 offset = binaryBlockStart + chunkSizeOffeset(processors::globalSize());
-
-	for(auto i=0; i<procNo; i++)
-	{
-		offset += chunkSizes[i];
-	}
-
-	if( data.size() != toRecv )
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-
-	
-	MPI_Status status;
-	
-	if( MPI_File_read_at_all(
-		fh,
-		offset,
-		data.data(),
-		data.size(),
-		MPI_UNSIGNED_CHAR,
-		&status) != MPI_SUCCESS )
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-	
-	MPI_File_close(&fh);
-
-	uint64 lastPos = offset + data.size();
-	/*if( MPI_Allreduce(
-		&lastPos,
-		&lastPosRead_, 
-		1, 
-		MPI_UINT64_T,  
-		MPI_MAX, 
-		MPI_COMM_WORLD) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction<<
-		"Error in max_reduction for lastPosRead_"<<endl;
-		return false;
-	}*/
-
-	return maxReduction(lastPos, lastPosRead_);
-#else
-	return readDataSTD(wordPath, chunkSizes, data, binaryBlockStart);
-#endif
-
-	
-}
-
-
-bool pFlow::dataIO::readMetaMPI
-(
-	const word& wordPath, 
-	std::vector<uint64>& chunkSizes,
-	uint64 startPosSearch,
-	uint64 &startPosBinaryBlock
-)
-{
-
-#ifdef pFlow_Build_MPI	
-	MPI_File fh;
-
-	if(MPI_File_open(
-		MPI_COMM_WORLD, 
-		wordPath.c_str(),
-		MPI_MODE_RDONLY,
-		MPI_INFO_NULL,
-		&fh) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-	
-
-	uint64 startPos = findBindaryBlockFlagMPI(fh, startPosSearch);
-	if( startPos == pFlow::dataIO::ErrorReturn )
-	{
-		fatalErrorInFunction;
-		return false;	
-	} 
-
-	startPosBinaryBlock  = startPos; 
-
-	uint64 numProcInFile;
-
-	if( MPI_File_read_at_all(
-		fh,
-		startPos,
-		&numProcInFile,
-		1,
-		MPI_UINT64_T,
-		MPI_STATUS_IGNORE) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-
-	
-	chunkSizes.resize(numProcInFile);
-	
-	if(MPI_File_read_at_all(
-		fh,
-		startPos + sizeof(numProcInFile),
-		chunkSizes.data(), 
-		chunkSizes.size(), 
-		MPI_UINT64_T,
-		MPI_STATUS_IGNORE) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-	MPI_File_close(&fh);
-
-	uint64 lastPos = startPos + sizeof(numProcInFile) + chunkSizes.size();
-	return maxReduction(lastPos, lastPosRead_);
-
-	/*if( MPI_Allreduce(
-		&lastPos,
-		&lastPosRead_, 
-		1, 
-		MPI_UINT64_T,  
-		MPI_MAX, 
-		MPI_COMM_WORLD) != MPI_SUCCESS)
-	{
-		fatalErrorInFunction<<
-		"Error in max_reduction for lastPosRead_"<<endl;
-		return false;
-	}*/
-
-#else
-	return readMetaSTD(wordPath, chunkSizes, startPosBinaryBlock);
-#endif
-}
-
-bool pFlow::dataIO::readMetaSTD
-(
-	const word& wordPath, 
-	std::vector<uint64>& chunkSizes,
-	uint64 startPosSearch,
-	uint64 &startPosBinaryBlock
-)
-{
-
-	std::FILE *fh = std::fopen(wordPath.c_str(), "rb");
-
-	if(!fh)
-	{
-		fatalErrorInFunction<<
-		"Error in Opening file "<< wordPath<<endl;
-		return false;
-	}
-
-	// set the start position for search 
-	std::fseek(fh, startPosSearch, SEEK_SET); 
-
-
-	uint64 startPos = findBindaryBlockFlagSTD(fh);
-	if(startPos == ErrorReturn )
-	{
-		fatalErrorInFunction;
-		return false;
-	}
-
-	startPosBinaryBlock = startPos;
-
-
-	std::fseek(fh, startPos, SEEK_SET); 
-
-	uint64 numProcInFile;
-	
-	auto res = std::fread(&numProcInFile, sizeof(numProcInFile), 1, fh);
-	
-
-	if(res != 1 )
-	{
-		fatalErrorInFunction<<
-		"Error in reading file "<< wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-
-	
-	chunkSizes.resize(numProcInFile);
-
-	res = std::fread(chunkSizes.data(), sizeof(numProcInFile), numProcInFile, fh);
-
-	if(res!= numProcInFile)
-	{
-		fatalErrorInFunction<<
-		"Error in reading chunkSizes from file "<< wordPath<<endl;
-		std::fclose(fh);
-		return false;
-	}
-
-	lastPosRead_ = static_cast<uint64>(std::ftell(fh));
-	
-	std::fclose(fh);
-	return true;
-}
-
-bool pFlow::dataIO::waitForAllMPI()
-{
-#ifdef pFlow_Build_MPI
-	MPI_Barrier(MPI_COMM_WORLD);
-#endif
-
-	return true;
-}
-
-
-bool pFlow::dataIO::maxReduction( uint64& src, uint64& dst)
-{
-	
-#ifdef pFlow_Build_MPI
-	if(processors::isParallel())
-	{
-		if(MPI_Allreduce(
-			&src,
-			&dst,
-			1,
-			MPI_UINT64_T,
-			MPI_MAX,
-			MPI_COMM_WORLD) != MPI_SUCCESS )
-		{
-			fatalErrorInFunction<<
-			"Error in max_reduction for write operation"<<endl;
-			return false;
-		}
-		return true;
 	}
 	else
 	{
-		dst = src;
+		printKeys
-		return true;
+		( 
+			fatalError << "Ctor Selector "<< Yellow_Text(ioMode) << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		fatalExit;
 	}
-#else
+	return nullptr;
-	dst = src;
-	return true;
-#endif
-	
 }
-
-bool pFlow::dataIO::BcastPos(uint64 & pos)
-{
-#ifdef pFlow_Build_MPI
-	if(processors::isParallel())
-	{
-		if( MPI_Bcast(
-			& pos,
-			1,
-			MPI_UINT64_T,
-			processors::masterNo(),
-			MPI_COMM_WORLD)!=MPI_SUCCESS )
-		{
-			fatalErrorInFunction<<
-			"Error in Bcast position"<<endl;
-			return false;
-		}
-	}
-#endif
-	return true;
-}

src/phasicFlow/streams/dataIO/dataIO.hpp

@@ -7,6 +7,10 @@
 #include "types.hpp"
 #include "span.hpp"
 #include "IOPattern.hpp"
+#include "iOstream.hpp"
+#include "iIstream.hpp"
+#include "virtualConstructor.hpp"
+
 
 namespace pFlow
 {
@@ -24,14 +28,124 @@ public:
 	
 protected:
 
-	IOPattern ioPattern_;
+	IOPattern 	ioPattern_;
 
-	uint64 lastPosWrite_ = 0;
+	word 		executionMode_;
 
-	uint64 lastPosRead_  = 0;
+	span<char> 			buffer_;
 
+	template<typename T>
+	span<char> createSpan( span<T> sp )
+	{
+		return span<char>( reinterpret_cast<char*>(sp.data()), sp.size()*sizeof(T) );
+	}
+
+public:
+
+	/// Type info
+	TypeInfo("dataIO");
+
+	dataIO(IOPattern::IOType iotype, word exeMode)
+	:
+	ioPattern_(iotype),
+	executionMode_(exeMode)
+	{}
+
+
+	dataIO(const dataIO&) = default;
+
+	dataIO(dataIO &&) = default;
+
+	dataIO& operator = (const dataIO&) = default;
+
+	dataIO& operator = (dataIO&&) = default;
+
+	virtual ~dataIO() = default;
+
+	create_vCtor
+	(
+		dataIO,
+		word,
+		(
+			IOPattern::IOType iotype, 
+			word exeMode
+		),
+		(iotype, exeMode)
+	);
+
+	/// Gather data from all processors and put the results in 
+	/// buffer_
+	virtual bool gatherData(span<char> sp ) = 0;
+
+	/// Write data to the end of file from all processors.
+	/// This method should be called from all processors. 
+	template<typename T>
+	bool writeData(iOstream& os, span<T> data);
+
+	template<typename T>
+	bool readData(iIstream& is, std::vector<T>& data);
+
+	template<typename T>
+	bool readAscii(
+		iIstream& is,
+		std::vector<T>& vec	);
+
+	static
+	uniquePtr<dataIO> create(IOPattern::IOType iotype, word exeMode);
+
+
+};
 	
-	bool writeDataToFileEndSTD(
+
+template<>
+inline 
+bool dataIO::writeData<word>(iOstream& os, span<word> data)
+{
+	/*std::ostringstream ist();
+
+	for(auto i=0; i<data.size(); i++)
+	{
+		ist<< data[i]<<std::endl;
+	}
+
+	span<char> sp(ist.str().data(), ist.str().size());
+
+	if(!gatherData( sp) )
+	{
+		fatalErrorInFunction<<
+		"Error in gathering data for out stream "<< os.name()<<endl;
+		return false;
+	}
+
+	uint64 len = buffer_.size();
+	os<< len << endl;
+
+	if(!os.write(buffer_.data(), buffer_.size()))
+	{
+		fatalErrorInFunction<<
+    	"error in writing binary data to "<<os.name()<<endl;
+    	return false;
+	}*/
+
+	if( !data.writeASCII(os) )
+    {
+    	fatalErrorInFunction<<
+    	"error in writing ASCII data to "<<os.name()<<endl;
+    	return false;
+    }
+
+	return true;
+}
+
+}
+
+#include "dataIOTemplate.cpp"
+
+#endif
+
+
+
+/*bool writeDataToFileEndSTD(
 		const word& wordPath, 
 		const span<unsigned char>& data);
 
@@ -143,14 +257,4 @@ bool pFlow::dataIO::readData<pFlow::word>(
 		return false;
 	}
 	return true;
-}
+}*/
-
-
-}
-
-#include "dataIOTemplate.cpp"
-
-#endif
-
-
-

src/phasicFlow/streams/dataIO/dataIONoMPI.cpp

@@ -0,0 +1,20 @@
+
+#include "dataIONoMPI.hpp"
+
+
+pFlow::dataIONoMPI::dataIONoMPI
+(
+	IOPattern::IOType iotype,
+	word exeMode
+)
+:
+	dataIO(iotype, exeMode)
+{}
+
+
+bool pFlow::dataIONoMPI::gatherData( span<char> sp )
+{
+	/// in serial mode, no gathering is required. 
+	buffer_ = sp;
+	return true;
+}

src/phasicFlow/streams/dataIO/dataIONoMPI.hpp

@@ -0,0 +1,48 @@
+#ifndef __datIONoMPI_hpp__
+#define __datIONoMPI_hpp__
+
+#include "dataIO.hpp"
+
+namespace pFlow
+{
+
+
+class dataIONoMPI
+:
+	public dataIO
+{
+protected:
+
+public:
+
+	TypeInfo("dataIO<NoMPI>");
+
+	dataIONoMPI(IOPattern::IOType iotype, word exeMode);
+
+	dataIONoMPI(const dataIONoMPI&) = default;
+
+	dataIONoMPI(dataIONoMPI&&) = default;
+
+
+	dataIONoMPI& operator=(const dataIONoMPI&) = default;
+
+	dataIONoMPI& operator=(dataIONoMPI&&) = default;
+
+	~dataIONoMPI() = default;
+
+	add_vCtor
+	(
+		dataIO,
+		dataIONoMPI,
+		word
+	);
+
+	bool gatherData(span<char> sp ) override;
+	
+};
+
+
+}
+
+
+#endif

src/phasicFlow/streams/dataIO/dataIOTemplate.cpp

@@ -1,4 +1,175 @@
+
 template<typename T>
+bool pFlow::dataIO::writeData(iOstream& os, span<T> data)
+{
+	/// first gather data from all processors (if any)
+	if(!gatherData( createSpan(data) ) )
+	{
+		fatalErrorInFunction<<
+		"Error in gathering data for out stream "<< os.name()<<endl;
+		return false;
+	}
+
+	if( ioPattern_.thisProcWriteData())
+	{
+		if( os.isBinary() )
+		{
+			// first write the size of data 		    
+			uint64 len = buffer_.size()/sizeof(T);
+			os<< len << endl;
+
+			if(!os.write(buffer_.data(), buffer_.size()))
+			{
+				fatalErrorInFunction<<
+		    	"error in writing binary data to "<<os.name()<<endl;
+		    	return false;
+			}
+
+		}
+		else
+		{
+
+			/// cast the data into T
+			span<T> allData( 
+				reinterpret_cast<T*>(buffer_.data()), 
+				buffer_.size()/sizeof(T));
+		    
+		    if( !allData.writeASCII(os) )
+		    {
+		    	fatalErrorInFunction<<
+		    	"error in writing ASCII data to "<<os.name()<<endl;
+		    	return false;
+		    }
+		}
+
+		return os.check(FUNCTION_NAME);
+	}
+	else
+	{
+		return true;
+	}
+}
+
+template<typename T>
+bool pFlow::dataIO::readData
+(
+	iIstream& is, 
+	std::vector<T>& data
+)
+{
+	
+	data.clear();
+
+	if(ioPattern_.thisProcReadData())
+	{
+		if(is.isBinary())
+		{
+			// read length of data 			
+			token firstToken(is);
+			
+			size_t len = 0;
+			if( firstToken.isInt64())
+			{
+				len = firstToken.int64Token();
+				data.resize(len);	
+			}
+			else
+			{
+				fatalErrorInFunction<<
+				"expected length of vector in the stream "<<is.name()<<
+				"but found this "<< firstToken<<endl;
+				return false;
+			}
+
+			is.read(reinterpret_cast<char*>(data.data()), len*sizeof(T));	
+			
+			return is.check(FUNCTION_NAME);
+		}
+		else
+		{
+			
+			return readAscii(is, data);
+		}
+	}
+	else
+	{
+		return true;
+	}
+}
+
+template<typename T>
+bool pFlow::dataIO::readAscii
+(
+	iIstream& is,
+	std::vector<T>& vec
+)
+{
+
+	if( !ioPattern_.thisProcReadData() ) return true;
+	
+
+	is.fatalCheck(FUNCTION_NAME);
+
+	vec.clear();
+	
+	token firstToken(is);
+
+	size_t len = 0;
+	if( firstToken.isInt64())
+	{
+		len = firstToken.int64Token();
+		vec.reserve(len);
+		firstToken  = token(is);
+	}
+
+	T val{};
+	if( firstToken.isPunctuation() ) // start of vector 
+	{
+		if(firstToken != token::BEGIN_LIST)
+		{
+			warningInFunction
+			<< "expected token "<< token::BEGIN_LIST 
+			<< " but found "<< firstToken ;
+			return false;
+
+		}
+
+		token lastToken(is);
+		
+		is.fatalCheck(FUNCTION_NAME);
+		
+		while(!(lastToken.isPunctuation()
+             	&& lastToken == token::END_LIST ))
+		{
+
+			is.putBack(lastToken);
+			
+			is >> val;
+			vec.push_back(val);
+			
+			is >> lastToken;
+			is.fatalCheck(FUNCTION_NAME);
+		}
+
+	} else
+	{
+		warningInFunction
+		<< "expected token "<< token::BEGIN_LIST 
+		<< " but found "<< firstToken<<endl; ;
+		return false;
+	}
+	
+	if(len>0&& len != vec.size())
+	{
+		warningInFunction<<"vector lendth specified "<< len << 
+		" is different from number of elements "<< vec.size()<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+/*template<typename T>
 bool pFlow::dataIO::writeDataEnd(
 	const word& wordPath,
 	const span<T>& data)
@@ -314,7 +485,7 @@ bool pFlow::dataIO::writeData
 	}
 	return true; 
 
-}
+}*/
 
 
 

src/phasicFlow/streams/iStream/iOstream.hpp

@@ -36,6 +36,13 @@ const inline char* whiteColor   = "\033[37m";
 
 const inline char* boldChar     = "\033[1m";
 
+#define Red_Text(text)      redColor<<text<<defaultColor
+#define Yellow_Text(text)   yellowColor<<text<<defaultColor
+#define Blue_Text(text)     blueColor<<text<<defaultColor
+#define Green_Text(text)    greenColor<<text<<defaultColor
+#define Magenta_Text(text)  magentaColor<<text<<defaultColor
+#define Cyan_Text(text)         cyanColor<<text<<defaultColor
+#define Bold_Text(text)         boldChar<<text<<defaultColor
 
 namespace pFlow
 {

src/phasicFlow/streams/streams.hpp

@@ -32,13 +32,6 @@ namespace pFlow
 
 }
 
-#define Red_Text(text)     	redColor<<text<<defaultColor
-#define Yellow_Text(text)  	yellowColor<<text<<defaultColor
-#define Blue_Text(text)    	blueColor<<text<<defaultColor
-#define Green_Text(text)   	greenColor<<text<<defaultColor
-#define Magenta_Text(text) 	magentaColor<<text<<defaultColor
-#define Cyan_Text(text) 		cyanColor<<text<<defaultColor
-#define Bold_Text(text) 		boldChar<<text<<defaultColor
 
 #define INFORMATION pFlow::mOutput<<boldChar<<magentaColor<<"> INFO: "<<defaultColor<<magentaColor
 #define END_INFO defaultColor<<pFlow::endl

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp

@@ -0,0 +1,23 @@
+
+
+
+#include "boundaryBase.hpp"
+
+
+pFlow::boundaryBase::boundaryBase
+(
+	const plane& 		bplane,
+	uint32 				mirrorProc, 
+	const word& 		name,
+	const word& 		type,
+	internalPoints& 	internal
+)
+:
+boundaryPlane_(bplane),
+name_(name),
+type_(type),
+mirrorProcessoNo_(mirrorProc),
+internal_(internal)
+{
+	
+}

src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp

@@ -0,0 +1,74 @@
+
+#ifndef __boundaryBase_hpp__
+#define __boundaryBase_hpp__
+
+
+#include "VectorSingles.hpp"
+#include "plane.hpp"
+
+
+
+namespace pFlow
+{
+
+// forward 
+
+class internalPoints;
+
+class boundaryBase
+{
+public:
+
+	enum DIRECTION: int8
+	{
+		
+	};
+
+protected:
+
+	plane 		boundaryPlane_;
+
+	word 		name_;
+
+	word 		type_;
+
+	uint32 		mirrorProcessoNo_;
+
+	/// list of particles indices on device 
+	uint32Vector_D 	indexList_;
+
+	/// a reference to 
+	internalPoints& internal_;
+
+public:
+
+	TypeInfo("boundaryBase");
+
+	boundaryBase(
+		const plane& 		bplane,
+		uint32 				mirrorProc, 
+		const word& 		name,
+		const word& 		type,
+		internalPoints& 	internal);
+
+	boundaryBase(const boundaryBase&) = default;
+
+	boundaryBase& operator=(const boundaryBase&) = default;
+
+	boundaryBase(boundaryBase&&) = delete;
+
+	boundaryBase& operator=(boundaryBase&&) = delete;
+
+	virtual ~boundaryBase() = default;
+
+	virtual bool update() = 0 ;
+
+	virtual bool iterate(uint32 iterNum, real t) = 0;
+
+	
+};
+
+}
+
+
+#endif //__boundaryBase_hpp__

src/phasicFlow/structuredData/boundaries/boundaryBase/scatterFieldAccess.hpp

@@ -0,0 +1,89 @@
+
+#ifndef __scatterFieldAccess_hpp__
+#define __scatterFieldAccess_hpp__
+
+
+#include "phasicFlowKokkos.hpp"
+
+namespace pFlow
+{
+
+
+template<typename T, typename ExecutionSpace>
+class scatterFieldAccess
+{
+public:
+	
+	using execution_space 	= ExecutionSpace;
+
+	using memory_space 		= typename execution_space::memory_space;
+
+protected:
+
+	uint32 								size_ = 0;
+
+	ViewType1D<uint32, memory_space> 	indices_;
+
+	ViewType1D<T, memory_space> 		fieldVals_;
+
+public:
+
+	scatterFieldAccess(
+		uint32 	sz
+		ViewType1D<uint32, memory_space> ind,
+		ViewType1D<T, memory_space> fVals)
+	:
+		size_(sz),
+		indices_(ind),
+		fieldVals_(fVals)
+	{}
+
+	scatterFieldAccess(const scatterFieldAccess&) = default;
+
+	scatterFieldAccess(scatterFieldAccess&&) = default;
+
+	scatterFieldAccess& operator=(const scatterFieldAccess&) = default;
+
+	scatterFieldAccess& operator=(scatterFieldAccess&&) = default;
+
+	~scatterFieldAccess() = default;
+
+	// - Methods
+
+	T& operator()(uint32 i)
+	{
+		return fieldVals_(indices_(i));
+	}
+
+	const T& operator()(uint32 i)const
+	{
+		return fieldVals_(indices_(i));
+	}
+
+	T& operator[](uint32 i)
+	{
+		return fieldVals_(indices_(i));
+	}
+
+	const T& operator[](uint32 i)const
+	{
+		return fieldVals_(indices_(i));
+	}
+
+	uint32 	size()const
+	{
+		return size_;
+	}
+
+	bool empty()const
+	{
+		return size_ == 0;
+	}
+
+};
+
+}
+
+
+#endif //__scatterFieldAccess_hpp__
+

src/phasicFlow/structuredData/box/box.hpp

@@ -62,18 +62,19 @@ public:
 		FUNCTION_H
 		box(iIstream& is);
 
-		FUNCTION_HD
+		INLINE_FUNCTION_HD
 		box(const box&) = default;
 
-		FUNCTION_HD
+		INLINE_FUNCTION_HD
 		box(box&&) = default;
 
-		FUNCTION_HD
+		INLINE_FUNCTION_HD
 		box& operator=(const box&) = default;
 
-		FUNCTION_HD
+		INLINE_FUNCTION_HD
 		box& operator=(box&&) = default;
 
+		INLINE_FUNCTION_HD
 		~box()=default;
 
 	//// - Methods 
@@ -85,13 +86,13 @@ public:
 		}
 
 		INLINE_FUNCTION_HD
-		realx3 minPoint()const
+		const realx3& minPoint()const
 		{
 			return min_;
 		}
 
 		INLINE_FUNCTION_HD
-		realx3 maxPoint()const
+		const realx3& maxPoint()const
 		{
 			return max_;
 		}

src/phasicFlow/structuredData/domain/domain.cpp

@@ -0,0 +1,59 @@
+/*------------------------------- 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 "domain.hpp"
+
+FUNCTION_H 
+pFlow::domain::domain(const box& db)
+:
+	domainBox_(db),
+	left_(
+		db.minPoint(),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.minPoint().z()),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.maxPoint().z()),
+		realx3(db.minPoint().x(), db.minPoint().y(), db.maxPoint().z())),
+	right_(
+		db.maxPoint(),
+		realx3(db.maxPoint().x(), db.maxPoint().y(), db.minPoint().z()),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.minPoint().z()),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.maxPoint().z())),
+	bottom_(
+		db.minPoint(),
+		realx3(db.minPoint().x(), db.minPoint().y(), db.maxPoint().z()),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.maxPoint().z()),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.minPoint().z())),
+	top_(
+		db.maxPoint(),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.maxPoint().z()),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.minPoint().z()),
+		realx3(db.maxPoint().x(), db.maxPoint().y(), db.minPoint().z())),
+	rear_(
+		db.minPoint(),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.minPoint().z()),
+		realx3(db.maxPoint().x(), db.maxPoint().y(), db.minPoint().z()),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.minPoint().z())),
+	front_(
+		db.maxPoint(),
+		realx3(db.maxPoint().x(), db.minPoint().y(), db.maxPoint().z()),
+		realx3(db.minPoint().x(), db.minPoint().y(), db.maxPoint().z()),
+		realx3(db.minPoint().x(), db.maxPoint().y(), db.maxPoint().z()))
+{
+
+}

src/phasicFlow/structuredData/domain/domain.hpp

@@ -0,0 +1,163 @@
+/*------------------------------- 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 __domain_hpp__
+#define __domain_hpp__
+
+
+#include "box.hpp"
+#include "plane.hpp"
+
+namespace pFlow
+{
+
+class domain
+{
+protected:
+	
+	box 		domainBox_;
+
+	/// -x
+	plane 		left_;
+
+	/// +x
+	plane 		right_;
+
+	/// -y
+	plane 		bottom_;
+
+	/// +y
+	plane		top_;
+
+	/// -z
+	plane 		rear_;
+
+	/// +z
+	plane 		front_;
+
+
+public:
+
+	// - type info
+	TypeInfoNV("domain");
+
+	//// - Constructors 
+		INLINE_FUNCTION_HD 
+		domain(){}
+
+		FUNCTION_H 
+		domain(const box& db);
+		
+		INLINE_FUNCTION_HD
+		domain(const domain&) = default;
+
+		INLINE_FUNCTION_HD
+		domain(domain&&) = default;
+
+		INLINE_FUNCTION_HD
+		domain& operator=(const domain&) = default;
+
+		INLINE_FUNCTION_HD
+		domain& operator=(domain&&) = default;
+
+		INLINE_FUNCTION_HD
+		~domain()=default;
+
+	//// - Methods 
+
+		INLINE_FUNCTION_HD
+		const auto& domainBox()const
+		{
+			return domainBox_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& left()const
+		{
+			return left_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& right()const
+		{
+			return right_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& bottom()const
+		{
+			return bottom_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& top()const
+		{
+			return top_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& rear()const
+		{
+			return rear_;
+		}
+
+		INLINE_FUNCTION_HD
+		const auto& front()const
+		{
+			return front_;
+		}
+
+	
+}; // domain
+
+INLINE_FUNCTION_HD
+bool equal(const domain& d1, const domain& d2)
+{
+	return equal(d1.domainBox(), d2.domainBox());
+}
+
+INLINE_FUNCTION_HD
+bool equal(const domain& d1, const domain& d2, real tol)
+{
+	return equal(d1.domainBox(), d2.domainBox(), tol);
+}
+
+INLINE_FUNCTION_HD
+bool operator ==(const domain& d1, const domain& d2)
+{
+	return equal(d1, d2);
+}
+
+/*FUNCTION_H
+iIstream& operator >>(iIstream& is, box& b);
+
+FUNCTION_H
+iOstream& operator << (iOstream& os, const box& b);
+
+INLINE_FUNCTION_HD
+box extendBox(const box& b, const realx3& dl)
+{
+	return box(b.minPoint()-dl , b.maxPoint()+dl);
+}*/
+
+}
+
+
+#endif

src/phasicFlow/structuredData/infinitePlane/infinitePlane.hpp

@@ -103,6 +103,21 @@ public:
 			return pointFromPlane(p)>=0;
 		}
 
+		INLINE_FUNCTION_HD
+		bool inPositiveDistance(const realx3& p, real dist)const
+		{	
+			real d = pointFromPlane(p);
+			return d >= 0.0 && d <= dist;
+		}
+
+		INLINE_FUNCTION_HD
+		bool inNegativeDistance(const realx3& p, real dist)const
+		{	
+			real d = pointFromPlane(p);
+			return d < 0.0 && d <= -dist;
+		}
+
+
 		INLINE_FUNCTION_HD
 		bool pointInNegativeSide(const realx3& p)const
 		{

src/phasicFlow/structuredData/pointStructure/internalPoints.cpp

@@ -0,0 +1,250 @@
+/*------------------------------- 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 "internalPoints.hpp"
+
+
+void pFlow::internalPoints::syncPFlag()const
+{
+	if(!pFlagSync_)
+	{
+		pFlagsH_ = pFlagsD_.clone<DefaultHostExecutionSpace>();
+		pFlagSync_ = true;
+	}
+}
+
+/*#include "setFieldList.hpp"
+#include "error.hpp"
+#include "iOstream.hpp"
+
+#include "mortonIndexing.hpp"*/
+
+/*FUNCTION_H
+bool pFlow::internalPoints::evaluateinternalPoints()
+{
+	if(pointFlag_.size() != pointPosition_.size())
+	{
+		fatalErrorInFunction<<
+		"number of elements in pointFlag and pointPosition is not equal \n";
+		return false;
+	}
+
+	setNumMaxPoints();
+
+	int32 minActive, maxActive;
+	numActivePoints_ = pFlow::internalPointsKernels::scanPointFlag(
+		0,
+		numPoints_,
+		static_cast<int8>(internalPoints::ACTIVE),
+		pointFlag_.deviceVectorAll(),
+		minActive,
+		maxActive
+		);
+	
+	activeRange_ = {minActive, maxActive};	
+	
+	return true;
+}
+
+FUNCTION_H
+void pFlow::internalPoints::setNumMaxPoints()
+{
+	maxPoints_ = pointFlag_.capacity();
+	numPoints_ = pointFlag_.size();
+}
+
+FUNCTION_H
+pFlow::realx3Field_D& pFlow::internalPoints::pointPosition()
+{
+	return pointPosition_;
+}
+
+FUNCTION_H
+pFlow::int8Field_HD& pFlow::internalPoints::pointFlag()
+{
+	return pointFlag_;
+}
+		
+pFlow::uniquePtr<pFlow::int32IndexContainer> 
+	pFlow::internalPoints::getNewPointsIndices(int32 numNewPoints)const
+{
+
+
+	if( capacity() - activeRange_.second >=  numNewPoints )
+	{
+		// fill the sequence starting from activeRange_.second-1
+		return makeUnique<int32IndexContainer>(
+			activeRange_.second,
+			activeRange_.second+numNewPoints);
+
+	}
+
+	// second, check if there is space at the beginning
+	if( activeRange_.first >= numNewPoints)
+	{
+		return makeUnique<int32IndexContainer>(
+			0,
+			numNewPoints);	
+	}
+
+	// otherwise scan the points from first to the end to find empty spaces 
+	int32Vector newPoints(
+		numNewPoints,
+		RESERVE());
+
+	newPoints.clear();
+	int32 numAdded = 0;
+	ForAll(i, pointFlag_)
+	{
+		if(!isActive(i))
+		{
+			newPoints.push_back(static_cast<int32>(i));
+			numAdded++;
+		}
+
+		if(numAdded == numNewPoints)
+		{
+			return makeUnique<int32IndexContainer>(
+				newPoints.data(),
+				numNewPoints);
+		}
+	}
+
+	// check if there is space at the end for the remaining of points
+	if( numAdded <numNewPoints && capacity() - size() >= numNewPoints - numAdded )
+	{
+		int32 ind = size();
+		for(int32 i=numAdded; i<numNewPoints; i++)
+		{
+			newPoints.push_back(ind);
+			ind++;
+		}
+		
+		return makeUnique<int32IndexContainer>(
+			newPoints.data(),
+			numNewPoints);
+	}
+	else
+	{
+		fatalErrorInFunction<<"not enough capacity for inserting particles into the point structure\n";
+		return nullptr;
+	}
+
+	return nullptr;
+}*/
+
+pFlow::internalPoints::internalPoints()
+:
+	pointPosition_("internalPoints", "internalPoints", initialCapacity_, 0, RESERVE()),
+	pFlagsD_(initialCapacity_, 0 , 0)
+{
+	syncPFlag();
+}
+
+
+pFlow::internalPoints::internalPoints
+(
+	const std::vector<realx3>& posVec
+)
+:
+	pointPosition_("internalPoints", "internalPoints", posVec.capacity(), 0, RESERVE()),
+	pFlagsD_(posVec.capacity(), 0, posVec.size())
+{
+
+	pointPosition_.assign(posVec);
+
+	syncPFlag();
+
+}
+
+
+const pFlow::internalPoints::pFlagTypeDevice& 
+	pFlow::internalPoints::activePointsMaskD() const 
+{
+	return pFlagsD_;
+}
+
+
+const pFlow::internalPoints::pFlagTypeHost&
+	pFlow::internalPoints::activePointsMaskH() const
+{
+	syncPFlag();
+	return pFlagsH_;
+}
+
+FUNCTION_H
+void pFlow::internalPoints::updateFlag
+(
+	const domain& dm, 
+	real dist
+)
+{
+	pFlagsD_.markDeleteInDomain
+	(
+		dm,
+		pointPosition_.deviceVectorAll(),
+		dist
+	);
+}
+
+
+FUNCTION_H
+bool pFlow::internalPoints::read
+(
+	iIstream& is, 
+	IOPattern::IOType iotype
+)
+{
+	
+	if( !pointPosition_.read(is, iotype))
+	{
+		fatalErrorInFunction<<
+		"Error in reading pointPosition from stream "<< is.name()<<endl;
+		return false;
+	}
+
+	pFlagsD_ = pFlagTypeDevice(pointPosition_.capacity(), 0, pointPosition_.size());
+	pFlagSync_ = false;
+	syncPFlag();
+
+	return true;
+}
+		
+
+
+FUNCTION_H
+bool pFlow::internalPoints::write
+(
+	iOstream& os, 
+	IOPattern::IOType iotype
+)const
+{
+	if( pFlagsD_.isAllActive())
+	{
+		return pointPosition_.write(os, iotype);
+	}
+	else
+	{
+		return pointPosition_.write(os, iotype, activePointsMaskH());
+	}
+}
+
+

src/phasicFlow/structuredData/pointStructure/internalPoints.hpp

@@ -0,0 +1,332 @@
+/*------------------------------- 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 __internalPoints_hpp__
+#define __internalPoints_hpp__
+
+#include "Fields.hpp"
+#include "pointFlag.hpp"
+
+
+//#include "indexContainer.hpp"
+
+
+
+namespace pFlow
+{
+
+//forward
+/*class box;
+class setFieldList;
+class repository;*/
+
+class internalPoints
+
+{
+public:
+	
+	inline static const uint32 initialCapacity_ = 10000;
+
+	using pointsType 		= realx3Field_D;
+
+	using device_type 		= typename pointsType::deviceType;
+
+	using memory_space 		= typename pointsType::memorySpace;
+
+	using execution_space 	= typename pointsType::executionSpace;
+
+	using pFlagTypeDevice 	= pointFlag<execution_space>;
+
+	using pFlagTypeHost 	= pointFlag<DefaultHostExecutionSpace>;
+
+
+protected:
+
+	//// - data members
+
+		// position of points on device
+		realx3Field_D 			pointPosition_;
+
+		// flag of points on device
+		mutable pFlagTypeDevice 		pFlagsD_;
+
+		mutable pFlagTypeHost 			pFlagsH_;
+
+		mutable bool 					pFlagSync_ = false;
+	
+
+	void syncPFlag()const;
+
+public:
+
+	friend class dynamicinternalPoints;
+	
+	// - type info
+	TypeInfo("internalPoints");
+	
+
+	//// - Constructors 
+
+		// - an empty internalPoints, good for reading from file 
+		internalPoints();
+
+		/// Construct from components
+		//internalPoints(const int8Vector& flgVec, const realx3Vector& posVec);
+
+		/// Construct from point positions, assume all points are active
+		internalPoints(const std::vector<realx3>& posVec);
+
+		/// Copy construct
+		internalPoints(const internalPoints&) = default;
+
+		/// Move construct 
+		internalPoints(internalPoints&&) = default;
+
+		/// Copy assignment 
+		internalPoints& operator=(const internalPoints&) = default;
+
+		/// Move assignment 
+		internalPoints& operator=(internalPoints&&) = delete;
+
+		/// Destructor 
+		virtual ~internalPoints() = default;
+
+
+	//// - Methods
+		FUNCTION_H
+		const pFlagTypeDevice& activePointsMaskD()const;
+
+
+		FUNCTION_H
+		const pFlagTypeHost& 	activePointsMaskH()const;
+
+		// - Const access pointPosition
+		FUNCTION_H
+		const realx3Field_D& pointPosition()const;
+		
+		
+
+		/*INLINE_FUNCTION_H
+		auto pointPositionHostAll()
+		{
+			return pointPosition_.hostVectorAll();
+		}*/
+
+		// - size of data structure
+		INLINE_FUNCTION_H
+		uint32 size()const
+		{
+			return pointPosition_.size();
+		}
+		
+		// - maximum capacity of data structure 
+		INLINE_FUNCTION_H
+		uint32 capacity()const
+		{
+			return pointPosition_.capacity();
+		}
+		
+		// - number of active points 
+		INLINE_FUNCTION_H
+		uint32 numActive() const
+		{
+			return pFlagsD_.numActive();
+		}
+				
+		// - if all points are active 
+		INLINE_FUNCTION_H
+		bool isAllActive()const
+		{
+			return pFlagsD_.isAllActive();
+		}
+
+		INLINE_FUNCTION_H
+		auto activeRange()const
+		{
+			return pFlagsD_.activeRange();
+		}
+
+		
+		FUNCTION_H
+		void updateFlag(const domain& dm, real dist);
+
+		/*FUNCTION_H
+		size_t markDeleteOutOfBox(const box& domain);*/
+
+		
+		///////////////////////////////////////////////////////////////////////////////////////////////////
+		
+		// - const access to points to be newly inserted
+		/*FUNCTION_H 
+		auto insertedPointIndex()const
+		{
+			return tobeInsertedIndex_;
+		}
+
+		FUNCTION_H 
+		auto insertedPointIndexH()const
+		{
+			return tobeInsertedIndex_.hostView();
+		}
+
+		FUNCTION_H 
+		auto insertedPointIndexD()const
+		{
+			return tobeInsertedIndex_.deviceView();
+		}
+
+
+		FUNCTION_H
+		auto mortonSortedIndex()const
+		{
+			return mortonSortedIndex_;
+		}
+
+
+		// - update data structure by inserting/setting new points 
+		//   Notifies all the fields in the registered list of data structure
+		//   and exclude the fields that re in the exclusionList
+		//   retrun nullptr if it fails 
+		/*FUNCTION_H
+		virtual uniquePtr<int32IndexContainer> insertPoints(
+			const realx3Vector& pos,
+			const setFieldList& setField,
+			repository& owner,
+			const List<eventObserver*>& exclusionList={nullptr}
+		);*/
+
+	
+	//// - IO operations 
+		
+		/// Read 
+		FUNCTION_H
+		bool read(iIstream& is, IOPattern::IOType iotype);
+		
+
+		/// Write 
+		FUNCTION_H
+		bool write(iOstream& os, IOPattern::IOType iotype)const;	
+
+};
+
+iOstream& operator<<(iOstream& os, const internalPoints& ip)
+{
+	if( !ip.write(os, IOPattern::AllProcessorsDifferent) )
+	{
+		ioErrorInFile(os.name(), os.lineNumber());
+		fatalExit;
+	}
+
+	return os;
+}
+
+} // pFlow
+
+
+
+
+#endif //__internalPoints_hpp__
+
+
+/*class activePointsDevice
+	{
+	protected:
+		ViewType1D<int8> 	flag_;
+
+		bool allActive_;
+
+		range activeRange_;
+
+	public:
+
+		INLINE_FUNCTION_H
+		activePointsDevice(bool allActive, range active, const ViewType1D<int8>& flag)
+		:
+			flag_(flag),
+			allActive_(allActive),
+			activeRange_(active)
+		{}
+
+		INLINE_FUNCTION_HD
+		activePointsDevice(const activePointsDevice&) = default;
+
+		INLINE_FUNCTION_HD
+		activePointsDevice& operator=(const activePointsDevice&) = default;
+
+		INLINE_FUNCTION_HD
+		bool operator()(int32 i)const {
+			if(i<activeRange_.second && flag_[i] == 1)return true;
+			return false;
+		}
+
+		INLINE_FUNCTION_HD
+		auto activeRange()const {
+			return activeRange_; 
+		}
+
+		INLINE_FUNCTION_HD
+		auto allActive()const {
+			return allActive_;
+		}
+
+	};
+
+	class activePointsHost
+	{
+	protected:
+		
+		ViewType1D<int8, HostSpace> 	flag_;
+
+		bool allActive_;
+
+		range activeRange_;
+
+	public:
+
+		INLINE_FUNCTION_H
+		activePointsHost(bool allActive, range active, const ViewType1D<int8, HostSpace>& flag)
+		:
+			flag_(flag),
+			allActive_(allActive),
+			activeRange_(active){}
+
+		INLINE_FUNCTION_H
+		activePointsHost(const activePointsHost&) = default;
+
+		INLINE_FUNCTION_H
+		activePointsHost& operator=(const activePointsHost&) = default;
+
+		INLINE_FUNCTION_H
+		bool operator()(int32 i)const {
+			if(i <activeRange_.second && flag_[i] == PointFlag::ACTIVE)return true;
+			return false;
+		}
+
+		INLINE_FUNCTION_H
+		auto activeRange()const{
+			return activeRange_;
+		}
+
+		INLINE_FUNCTION_H
+		bool allActive()const {
+			return allActive_;
+		}
+	};*/

src/phasicFlow/structuredData/pointStructure/pointFlag.hpp

@@ -0,0 +1,279 @@
+/*------------------------------- 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 __pointFlag_hpp__
+#define __pointFlag_hpp__
+
+
+#include "phasicFlowKokkos.hpp"
+#include "domain.hpp"
+
+
+
+namespace pFlow
+{
+
+
+template<typename ExecutionSpace>
+class pointFlag
+{
+	enum Flag: uint8
+	{
+		DELETED 	= 1, 	//10000000
+		INTERNAL 	= 2,	//00000001
+		LEFT 		= 4,	//00000010
+		RIGHT 		= 8,    //00000100
+		BOTTOM 		= 16,	//00001000
+		TOP 		= 32,	//00010000
+		REAR 		= 64, 	//00100000
+		FRONT 		= 128	//01000000
+	};
+
+	using execution_space 	= ExecutionSpace;
+
+	using memory_space 		= typename execution_space::memory_space;
+
+	using viewType 			= ViewType1D<uint8, memory_space>;
+
+	using device_type 		= typename viewType::device_type;
+
+protected:
+	
+	viewType 	flags_;
+
+	uint32 		numActive_ 		= 0;
+
+	rangeU32 	activeRange_ 	= {0,0};
+
+	bool 			isAllActive_	= false;
+
+	uint32 		nLeft_ 	= 0;
+
+	uint32 		nRight_ = 0;
+
+	uint32 		nBottom_= 0;
+
+	uint32 		nTop_ 	= 0;
+
+	uint32 		nRear_ 	= 0;
+
+	uint32 		nFront_ = 0;
+
+	//- Protected methods 
+
+	
+
+public:
+
+	friend class internalPoints;
+
+	pointFlag()
+	{}
+
+	pointFlag(uint32 capacity, uint32 start, uint32 end )
+	:
+		flags_("pointFlag", capacity ),
+		numActive_(end-start),
+		activeRange_(start, end),
+		isAllActive_(true)
+	{
+		fill(flags_, 0, capacity, static_cast<uint8>(Flag::DELETED));
+		fill(flags_, activeRange_, static_cast<uint8>(Flag::INTERNAL));
+	}
+
+	pointFlag(
+		viewType 	flags,
+		uint32 		numActive,
+		rangeU32 	activeRange,
+		bool 		isAllActive)
+	:
+		flags_(flags),
+		numActive_(numActive),
+		activeRange_(activeRange),
+		isAllActive_(isAllActive)
+	{}
+
+	pointFlag(const pointFlag&) = default;
+
+	pointFlag& operator=(const pointFlag&) = default;
+
+	pointFlag(pointFlag&&) = default;
+
+	pointFlag& operator=(pointFlag&&) = default;
+
+	~pointFlag() = default;	
+
+	INLINE_FUNCTION_HD
+	bool isAllActive()const
+	{
+		return isAllActive_;
+	}
+
+	INLINE_FUNCTION_HD
+	const auto& activeRange()const
+	{
+		return activeRange_;
+	}
+
+	INLINE_FUNCTION_HD
+	auto numActive()const
+	{
+		return numActive_;
+	}
+
+	INLINE_FUNCTION_HD
+	bool operator()(uint32 i)
+	{
+		return isActive(i);
+	}
+
+	INLINE_FUNCTION_HD
+	bool isActive(uint32 i)const
+	{
+		return flags_[i] > DELETED;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isBoundary(uint32 i)const
+	{
+		return flags_[i] > INTERNAL;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isBoundary(uint8 flg)const
+	{
+		return flg > INTERNAL;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isLeft(uint32 i)const
+	{
+		return (flags_[i]&LEFT) == LEFT;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isLeft(uint8 flg)const
+	{
+		return (flg&LEFT) == LEFT;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isRight(uint32 i)const
+	{
+		return (flags_[i]&&RIGHT) == RIGHT;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isRight(uint8 flg)const
+	{
+		return (flg&&RIGHT) == RIGHT;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isBottom(uint32 i)const
+	{
+		return (flags_[i]&&BOTTOM) == BOTTOM;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isBottom(uint8 flg)const
+	{
+		return (flg&&BOTTOM) == BOTTOM;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isTop(uint32 i)const
+	{
+		return (flags_[i]&&TOP) == TOP;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isTop(uint8 flg)const
+	{
+		return (flg&&TOP) == TOP;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isRear(uint32 i)const
+	{
+		return (flags_[i]&&REAR) == REAR;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isRear(uint8 flg)const
+	{
+		return (flg&&REAR) == REAR;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isFront(uint32 i)const
+	{
+		return (flags_[i]&&FRONT) == FRONT;
+	}
+
+	INLINE_FUNCTION_HD
+	bool isFront(uint8 flg)const
+	{
+		return (flg&&FRONT) == FRONT;
+	}
+
+
+	template<typename ExeSpace>
+	pointFlag<ExeSpace> clone()const
+	{
+		using newViewType = ViewType1D<uint8, typename ExeSpace::memory_space>;
+		newViewType newFlags(
+			"pointFlag",
+			flags_.size());
+
+		copy(newFlags, flags_);
+
+		return pointFlag<ExeSpace>(
+			newFlags,
+			numActive_,
+			activeRange_,
+			isAllActive_);
+	}
+
+	
+	uint32 scanPointFlag();
+
+	uint32 markDeleteInDomain(
+		domain 								dm,
+		ViewType1D<realx3, memory_space> 	points,
+		real dist);
+
+	void fillNeighborsLists(
+		ViewType1D<uint32, memory_space> leftList,
+		ViewType1D<uint32, memory_space> rightList,
+		ViewType1D<uint32, memory_space> bottomList,
+		ViewType1D<uint32, memory_space> topList,
+		ViewType1D<uint32, memory_space> rearList,
+		ViewType1D<uint32, memory_space> frontList);
+
+
+
+};
+
+
+}
+
+#include "pointFlagKernels.hpp"
+
+#endif // __pointFlag_hpp__

src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp

@@ -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.
+
+-----------------------------------------------------------------------------*/
+#ifndef __pointFlagKernels_hpp__
+#define __pointFlagKernels_hpp__
+
+
+
+template<typename ExecutionSpace>
+pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::scanPointFlag()
+{
+		
+	using rpScanFlag = Kokkos::RangePolicy<execution_space,  
+			Kokkos::IndexType<uint32>>;
+
+	uint32 numActive = 0;
+	
+	
+	uint32 start = activeRange().start();
+	uint32 end   = activeRange().end();
+
+	uint32 minRange = end;
+	uint32 maxRange = start;
+
+	if(start<end)
+	{
+		Kokkos::parallel_reduce(
+			"pointFlagKernels::scanPointFlag",
+			rpScanFlag(start, end),
+			CLASS_LAMBDA_HD(
+				uint32 i, 
+				uint32& minUpdate, 
+				uint32& maxUpdate, 
+				uint32& sumToUpdate)
+			{
+				if(this->isActive(i))
+				{
+					sumToUpdate++;
+					minUpdate = min(minUpdate,i);
+					maxUpdate = max(maxUpdate,i);
+				}
+			},
+			Kokkos::Min<uint32>(minRange), 
+			Kokkos::Max<uint32>(maxRange),
+			numActive);	
+	}
+	
+	if(numActive==0)
+	{
+		minRange = 0;
+		maxRange = 0;
+	}
+	else
+	{
+		// add one to maxRange to make it half-open
+		maxRange ++;
+	}
+
+	activeRange_ = {minRange, maxRange};
+	numActive_ 	= numActive;
+	isAllActive_ = activeRange_.numElements() == numActive_;
+
+	return numActive;
+}
+
+
+template<typename ExecutionSpace>
+pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markDeleteInDomain
+(
+	domain 								dm,
+	ViewType1D<realx3, memory_space> 	points,
+	real dist)
+{
+	
+	using rpMark = Kokkos::RangePolicy<execution_space,  
+			Kokkos::IndexType<uint32>>;
+
+	uint32 start = activeRange().start();
+	uint32 end   = activeRange().end();
+
+	uint32 minRange = end;
+	uint32 maxRange = start;
+
+	uint32 numMarked = 0;
+	uint32 nLeft = 0, nRight = 0;
+	uint32 nBottom = 0, nTop = 0;
+	uint32 nRear = 0, nFront = 0;
+
+	if(start<end)
+	{
+		
+		Kokkos::parallel_reduce(
+			"pointFlagKernels::markDeleteInDomain",
+			rpMark(start,end),
+			CLASS_LAMBDA_HD(
+				uint32 i, 
+				uint32& minUpdate, 
+				uint32& maxUpdate, 
+				uint32& valDelete,
+				uint32& valLeft,
+				uint32& valRight,
+				uint32& valBottom,
+				uint32& valTop,
+				uint32& valRear,
+				uint32& valFront){
+				if(this->isActive(i))
+				{
+					realx3 p = points[i];
+					if( !dm.domainBox().isInside(p) )
+					{
+						flags_[i] = DELETED;
+						valDelete++;
+					}
+					else
+					{
+
+						uint8 flg  = 0;
+						minUpdate = min(minUpdate,i);
+						maxUpdate = max(maxUpdate,i);
+
+						if(dm.left().inPositiveDistance(p, dist))
+						{
+							flg += LEFT;
+							valLeft++;	
+						} 
+
+						if(dm.right().inPositiveDistance(p, dist))
+						{
+							flg += RIGHT;
+							valRight++;
+						}
+						
+						if(dm.bottom().inPositiveDistance(p, dist))
+						{
+							flg += BOTTOM;
+							valBottom++;
+						} 
+
+						if(dm.top().inPositiveDistance(p, dist)) 
+						{
+							flg += TOP;
+							valTop++;
+						}
+						
+						if(dm.rear().inPositiveDistance(p, dist)) 
+						{
+							flg += REAR;
+							valRear++;
+						}
+
+						if(dm.front().inPositiveDistance(p, dist))
+						{
+							flg += FRONT;
+							valFront++;
+						} 
+
+						flags_[i] = flg>static_cast<uint8>(0)? flg: INTERNAL;
+					}
+				}
+			},
+			Kokkos::Min<uint32>(minRange), 
+			Kokkos::Max<uint32>(maxRange),
+			numMarked,
+			nLeft,
+			nRight,
+			nBottom,
+			nTop,
+			nRear,
+			nFront);
+
+	}
+	
+	// means either range was empty or all points have been deleted.
+	if(minRange<start || maxRange>end)
+	{
+		minRange = 0;
+		maxRange = 0;
+	}
+	else
+	{
+		maxRange++; // add one to make it half
+	}
+
+	activeRange_ = {minRange, maxRange};
+	isAllActive_ = isAllActive_ && numMarked == 0; 
+	numActive_ 	-= numMarked;
+
+	nLeft_ 	= nLeft;
+	nRight_	= nRight;
+	nBottom_= nBottom;
+	nTop_ 	= nTop;
+	nRear_	= nRear;
+	nFront_ = nFront;
+
+	return numMarked;
+}
+
+template<typename ExecutionSpace>
+void pFlow::pointFlag<ExecutionSpace>::fillNeighborsLists
+(
+	ViewType1D<uint32, memory_space> leftList,
+	ViewType1D<uint32, memory_space> rightList,
+	ViewType1D<uint32, memory_space> bottomList,
+	ViewType1D<uint32, memory_space> topList,
+	ViewType1D<uint32, memory_space> rearList,
+	ViewType1D<uint32, memory_space> frontList
+)
+{
+	using rpMark = Kokkos::RangePolicy<execution_space,  
+			Kokkos::IndexType<uint32>>;
+
+	uint32 start = activeRange().start();
+	uint32 end   = activeRange().end();
+
+	ViewType1D<uint32, memory_space> nElems("nElems",6);
+
+	fill(nElems, 0, 6, 0);
+
+	if(start<end)
+	{
+		Kokkos::parallel_for(
+			"pointFlagKernels::markDeleteInDomain",
+			rpMark(start,end),
+			CLASS_LAMBDA_HD(uint32 i){
+
+				uint32 flg = flags_(i);
+
+				if(this->isBoundary(flg))
+				{
+					if(this->isLeft(flg)) leftList[Kokkos::atomic_fetch_add(&nElems[0],1)] = i;
+					if(this->isRight(flg)) rightList[Kokkos::atomic_fetch_add(&nElems[1],1)] = i;
+					if(this->isBottom(flg)) bottomList[Kokkos::atomic_fetch_add(&nElems[2],1)] = i;
+					if(this->isTop(flg)) topList[Kokkos::atomic_fetch_add(&nElems[3],1)] = i;
+					if(this->isRear(flg)) rearList[Kokkos::atomic_fetch_add(&nElems[4],1)] = i;
+					if(this->isFront(flg)) frontList[Kokkos::atomic_fetch_add(&nElems[5],1)] = i;
+				}
+			});
+		Kokkos::fence();	
+	}
+	
+
+}
+
+#endif // __pointFlagKernels_hpp__

src/phasicFlow/structuredData/pointStructure/pointStructure.hpp

@@ -23,6 +23,8 @@ Licence:
 #define __pointStructure_hpp__
 
 
+#include "Lists.hpp"
+#include "internalPoints.hpp"
 
 #include "Vectors.hpp"
 #include "VectorSingles.hpp"
@@ -37,9 +39,8 @@ namespace pFlow
 {
 
 //forward
-class box;
+class boundaryBase;
-class setFieldList;
+
-class repository;
 
 class pointStructure
 :
@@ -47,152 +48,24 @@ class pointStructure
 {
 public:
 	
-	enum PointFlag: int8
+	using boundaryListType = ListPtr<boundaryBase>;
-	{
-		DELETED 	= -1,
-		ACTIVE 		= 1
-	};
-
-	
-	inline static const size_t maxSizeDefault_ = 10000;
-
-	class activePointsDevice
-	{
-	protected:
-		ViewType1D<int8> 	flag_;
-
-		bool allActive_;
-
-		range activeRange_;
-
-	public:
-
-		INLINE_FUNCTION_H
-		activePointsDevice(bool allActive, range active, const ViewType1D<int8>& flag)
-		:
-			flag_(flag),
-			allActive_(allActive),
-			activeRange_(active)
-		{}
-
-		INLINE_FUNCTION_HD
-		activePointsDevice(const activePointsDevice&) = default;
-
-		INLINE_FUNCTION_HD
-		activePointsDevice& operator=(const activePointsDevice&) = default;
-
-		INLINE_FUNCTION_HD
-		bool operator()(int32 i)const {
-			if(i<activeRange_.second && flag_[i] == 1)return true;
-			return false;
-		}
-
-		INLINE_FUNCTION_HD
-		auto activeRange()const {
-			return activeRange_; 
-		}
-
-		INLINE_FUNCTION_HD
-		auto allActive()const {
-			return allActive_;
-		}
-
-	};
-
-	class activePointsHost
-	{
-	protected:
-		
-		ViewType1D<int8, HostSpace> 	flag_;
-
-		bool allActive_;
-
-		range activeRange_;
-
-	public:
-
-		INLINE_FUNCTION_H
-		activePointsHost(bool allActive, range active, const ViewType1D<int8, HostSpace>& flag)
-		:
-			flag_(flag),
-			allActive_(allActive),
-			activeRange_(active){}
-
-		INLINE_FUNCTION_H
-		activePointsHost(const activePointsHost&) = default;
-
-		INLINE_FUNCTION_H
-		activePointsHost& operator=(const activePointsHost&) = default;
-
-		INLINE_FUNCTION_H
-		bool operator()(int32 i)const {
-			if(i <activeRange_.second && flag_[i] == PointFlag::ACTIVE)return true;
-			return false;
-		}
-
-		INLINE_FUNCTION_H
-		auto activeRange()const{
-			return activeRange_;
-		}
-
-		INLINE_FUNCTION_H
-		bool allActive()const {
-			return allActive_;
-		}
-	};
 
 protected:
 
 	//// - data members 
+		domain 				globalDomain_;
 
-		// number of points / size of structure
+		domain 				thisDomain_;
-		size_t 			numPoints_ 		= 0;	
 
-		// maximum number of points
+		internalPoints 		internal_;
-		size_t 			maxPoints_ 		= maxSizeDefault_;
 
-		// flag of points on device
+		boundaryListType 	boundaries_;
-		int8Field_HD 			pointFlag_;
 
-		// position of points on device
+		wordList 			boundaryTypes_;	
-		realx3Field_D 			pointPosition_;
-
-		// number of active points
-		size_t 					numActivePoints_ = 0;
-
-		// index range of active points (half-open range)
-		range 					activeRange_;	
-
-		/// Index vector for points to be inserted 
-		int32IndexContainer   	tobeInsertedIndex_;
-
-		/// Sorted index of particles based on morton code 
-		int32IndexContainer 	mortonSortedIndex_;
-	
-
-	//// - protected methods
-		FUNCTION_H
-		bool evaluatePointStructure();
-
-		FUNCTION_H
-		void setNumMaxPoints();
-
-		// - access to pointPosition
-		FUNCTION_H
-		realx3Field_D& pointPosition();
-		
-		// - access to pointFlag
-		FUNCTION_H
-		int8Field_HD& pointFlag();
-
-		FUNCTION_H
-		uniquePtr<int32IndexContainer> 
-			getNewPointsIndices(int32 numNewPoints)const;
 
 public:
 
-	friend class dynamicPointStructure;
+		
-	
 	// - type info
 	TypeInfo("pointStructure");
 	
@@ -228,41 +101,7 @@ public:
 		pointStructure& operator=(pointStructure&&) = delete;
 
 		// - destructor 
-		virtual ~pointStructure() = default;
+		virtual ~pointStructure() = default
-
-
-	//// - Methods
-		activePointsDevice activePointsMaskD()const
-		{
-			return activePointsDevice(
-				this->allActive(),
-				activeRange(),
-				pointFlag_.deviceVectorAll()
-				);
-		}
-
-		activePointsHost activePointsMaskH()const
-		{
-			return activePointsHost(
-				this->allActive(),
-				activeRange(),
-				pointFlag_.hostVectorAll()
-				);
-		}
-
-		// - const access pointPosition
-		FUNCTION_H
-		const realx3Field_D& pointPosition()const;
-		
-		// - const access to pointFlag
-		FUNCTION_H
-		const int8Field_HD& pointFlag()const;
-
-		INLINE_FUNCTION_H
-		auto pointPositionHostAll()
-		{
-			return pointPosition_.hostVectorAll();
-		}
 
 		// - size of data structure
 		FUNCTION_H
@@ -272,103 +111,37 @@ public:
 		FUNCTION_H
 		label capacity()const;
 		
-		// - number of active points 
-		FUNCTION_H
-		label numActive() const;
-				
-		// - if all points are active 
-		FUNCTION_H
-		bool allActive()const;
-
-		INLINE_FUNCTION_H
-		range activeRange()const
-		{
-			return activeRange_;
-		}
-
-		INLINE_FUNCTION_H
-		bool isActive(label i)const
-		{
-			return pointFlag_[i] == ACTIVE;
-		}
-
-		FUNCTION_H
-		size_t markDeleteOutOfBox(const box& domain);
-
-		// - update data structure by removing the marked points 
-		//   Notifies all the fields that are built based on this data structure
-		//   and then apply removing to the pointPosition_ and pointFlag_
-		FUNCTION_H
-		virtual bool updateForDelete();
-
-
-		FUNCTION_H
-		virtual bool mortonSortPoints(const box& domain, real dx);
-
-		///////////////////////////////////////////////////////////////////////////////////////////////////
 		
-		// - const access to points to be newly inserted
-		FUNCTION_H 
-		auto insertedPointIndex()const
-		{
-			return tobeInsertedIndex_;
-		}
-
-		FUNCTION_H 
-		auto insertedPointIndexH()const
-		{
-			return tobeInsertedIndex_.hostView();
-		}
-
-		FUNCTION_H 
-		auto insertedPointIndexD()const
-		{
-			return tobeInsertedIndex_.deviceView();
-		}
-
-
-		FUNCTION_H
-		auto mortonSortedIndex()const
-		{
-			return mortonSortedIndex_;
-		}
-
-
 		// - update data structure by inserting/setting new points 
 		//   Notifies all the fields in the registered list of data structure
 		//   and exclude the fields that re in the exclusionList
 		//   retrun nullptr if it fails 
 		FUNCTION_H
-		virtual uniquePtr<int32IndexContainer> insertPoints(
+		/*virtual uniquePtr<int32IndexContainer> insertPoints(
 			const realx3Vector& pos,
 			const setFieldList& setField,
 			repository& owner,
 			const List<eventObserver*>& exclusionList={nullptr}
-		);
+		);*/
+
+		boundaryBase& boundary(size_t i)
+		{
+			return boundaries_[i];
+		}
+
+
+		const boundaryBase& boundary(size_t i)const
+		{
+			return boundaries_[i];
+		}
+
 
 	
 	//// - IO operations 
-		// - read pointStructure from is
+		
-		FUNCTION_H
+		
-		bool readPointStructure(iIstream& is);
 
-		// - write pointStructure to os
-		FUNCTION_H
-		bool writePointStructure(iOstream& os)const;
 
-		// - read 
-		FUNCTION_H
-		bool read(iIstream& is)
-		{
-			return readPointStructure(is);
-		}
-
-		// - write 
-		FUNCTION_H
-		bool write(iOstream& os)const
-		{
-			return writePointStructure(os);
-		}
 
 };
 

src/phasicFlow/types/basicTypes/bTypesFunctions.hpp

@@ -190,9 +190,9 @@ bool readValue(const word& w, bool& val)
 
 
 INLINE_FUNCTION_HD 
-bool equal(const real& s1, const real& s2)
+bool equal(const real& s1, const real& s2, real tol = smallValue)
 {
-    return abs(s1 - s2) <= smallValue;
+    return abs(s1 - s2) <= tol;
 }
 
 INLINE_FUNCTION_HD 

src/phasicFlow/types/triple/triple.hpp

@@ -231,7 +231,17 @@ struct triple
 };
 
 template<typename T>
-bool INLINE_FUNCTION_HD equal( const triple<T>& opr1, const triple<T>& opr2 );
+bool INLINE_FUNCTION_HD equal( const triple<T>& opr1, const triple<T>& opr2);
+
+
+bool
+INLINE_FUNCTION_HD
+equal(const triple<real>& opr1, const triple<real>& opr2, real tol)
+{
+	return 	equal( opr1.x(), opr2.x(), tol ) && 
+			equal( opr1.y(), opr2.y(), tol ) && 
+			equal( opr1.z(), opr2.z(), tol );
+}
 
 
 } /// end of pFlow