Merge pull request #103 from PhasicFlow/develop

Particle insertion upgrade and pointStructure selectors

src/Interaction/contactSearch/methods/particleWallContactSearchs.cpp

@@ -34,7 +34,7 @@ pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
     domainBox_(domain),
     updateInterval_
     (
-        dict.getValOrSet<uint32>("updateInterval", 1)
+        max(dict.getValOrSet<uint32>("updateInterval", 1),1u)
     )
 {
     

src/Particles/Insertion/Insertion/Insertion.cpp

@@ -18,31 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-/*template<typename ShapeType>
-bool pFlow::Insertion<ShapeType>::readInsertionDict
-(
-	const dictionary& dict
-)
-{
-	if(!insertion::readInsertionDict(dict)) return false;
-	
-	regions_.clear();
-	
-	if( !this->isActive() )
-	{
-		return true;
-	}
 
-	wordList regionDicNames = dict.dictionaryKeywords();
-
-	for(auto& name:regionDicNames)
-	{
-		REPORT(2)<<"reading insertion region "<< greenText(name)<<endREPORT;
-		regions_.push_backSafe(dict.subDict(name), shapes_);
-	}
-
-	return true;
-}
 
 template<typename ShapeType>
 bool pFlow::Insertion<ShapeType>::writeInsertionDict
@@ -50,7 +26,8 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
 	dictionary& dict
 )const
 {
-	if( !insertion::writeInsertionDict(dict) ) return false;
+	
+	if(!insertion::writeInsertionDict(dict))return false;
 
 	if( !this->isActive() ) return true;
 
@@ -60,16 +37,18 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
 
 		if( !regions_[i].write(rgnDict) )
 		{
+			fatalErrorInFunction<<
+			"Error in writing to dictionary "<<rgnDict.globalName()<<endl;
 			return false;
 		}
 	}
 
 	return true;
-}*/
+}
 
 template<typename ShapeType>
 bool
-pFlow::Insertion<ShapeType>::setInsertionRegions()
+pFlow::Insertion<ShapeType>::readInsertionDict()
 {
 	regions_.clear();
 	
@@ -102,7 +81,7 @@ pFlow::Insertion<ShapeType>::Insertion(
 	insertion(prtcl),
 	shapes_(shapes)
 {
-	if(!setInsertionRegions())
+	if(!readInsertionDict())
 	{
 		fatalErrorInFunction;
 		fatalExit;

src/Particles/Insertion/Insertion/Insertion.hpp

@@ -62,12 +62,11 @@ private:
 	// - insertion regions
 	ListPtr<InsertionRegion<ShapeType>> regions_;
 
+  bool readInsertionDict();
 
-	bool setInsertionRegions();
+protected:
 
-	// bool readInsertionDict(const dictionary& dict);
-
-	// bool writeInsertionDict(dictionary& dict)const;
+	bool writeInsertionDict(dictionary& dict)const override;
 
 public:
 

src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp

@@ -55,8 +55,8 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
 	uint32 iter, 
 	real t,
 	real dt,
-	wordVector& names,
-	realx3Vector& pos,
+	wordVector& 	names,
+	realx3Vector& 	positions,
 	bool& insertionOccured
 )
 {
@@ -65,65 +65,185 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
 	if(!insertionTime(iter, t, dt)) return true;
 
 	uint32 newNum = numberToBeInserted(iter, t, dt);
-	
 	if(newNum == 0) return true;
 
+	// get the internal box 
 	auto internalBox =  pStruct().internalDomainBox();
-
 	if( !(internalBox.minPoint() < internalBox.maxPoint()))
 	{
 		WARNING<<"Minimum point of internal point is not lower than "<<
 		"the maximum point \n"<<
 		"minimum point: "<< internalBox.minPoint()<<
 		"\nmaximum point:"<<internalBox.maxPoint()<<END_WARNING;
+
+		return false;
 	}
 
-	names.reserve(max(newNum,names.capacity()));
-	pos.reserve(max(newNum,pos.capacity()));
-	names.clear();
-	pos.clear();
-
-	realVector diams("diams", newNum, 0, RESERVE());
-
 	
-	size_t n = 0;
-	uint32 idx;
+
+	/// reserve enough space 
+	names.reserve(max(newNum,names.capacity()));
+	names.clear();
+	
+	positions.reserve(max(newNum,positions.capacity()));
+	positions.clear();
+
+	auto allPositions = realx3Vector("allPositions");
+	auto allDiameters = realVector("allDiameters");
+
 	auto& mix = mixture();
 	auto& pReg = pRegion();
+
+	real maxDiam = shapes_.maxBoundingSphere();
+	auto minP = pReg.minPoint() - maxDiam;
+	auto maxP = pReg.maxPoint() + maxDiam;
+	
+	if(Insertion().checkForCollision())
+	{
+		// check for collision with already inserted particles 	
+		// so, use selector to select particles in the simulation 
+		auto bDict = dictionary("boxInfo");
+		bDict.add("min", minP);
+		bDict.add("max", maxP);
+		
+		auto selector = pStructSelector::create(
+			"box", 
+			pStruct(), 
+			bDict);
+		
+		allPositions = selector().selectedPointPositions();
+		allDiameters = selectedFieldVals<real>
+		(
+			selector(), 
+			Insertion().diameterName()
+		);
+	}
+	
+	auto collCheck = collisionCheck( 
+		{minP, maxP}, 
+		maxDiam,
+		allPositions,
+		allDiameters);
+
+	uint32 numInserted = 0;
+
+	uint32 idx;
 	word name = mix.getNextShapeName();
 	shapes_.shapeNameToIndex(name, idx);
 	real d = shapes_.boundingDiameter(idx);
 
 	for(uint32 i=0; i< 100*newNum ; ++i)
 	{
-		
 		realx3 p = pReg.peek();
-
 		// check if point is inside internal box 
 		if(!internalBox.isInside(p))continue;
 
-		if( !checkForContact(pos, diams, p, d) )
+		if( collCheck.checkPoint( p, d) )
 		{
 			names.push_back(name);
-			pos.push_back(p);
-			diams.push_back(d);
-			n++;
-
-			if( n == newNum ) 
-			{
-				addToNumInserted(newNum);
-				insertionOccured = true;
-				return true;
-			}
-
+			positions.push_back(p);
+			numInserted++;
+			if( numInserted == newNum ) break;
+			
+			// add this new particle to collision check set 
+			allDiameters.push_back(d);
+			allPositions.push_back(p);
+			collCheck.mapLastAddedParticle();
+							
+			// obtain next shape name and diameter 
 			name = mix.getNextShapeName();
 			shapes_.shapeNameToIndex(name, idx);
 			d = shapes_.boundingDiameter(idx);
 		}
-
+		
 	}
-
-	addToNumInserted(n);
+	
 	insertionOccured = true;
-	return false;
-}
+	addToNumInserted(numInserted);
+	return numInserted == newNum;
+}
+
+
+
+/*if(!checkForCollision)
+	{
+		realVector diams("diams", newNum, 0, RESERVE());
+		
+		uint32 idx;
+		word name = mix.getNextShapeName();
+		shapes_.shapeNameToIndex(name, idx);
+		real d = shapes_.boundingDiameter(idx);	
+
+		for(uint32 i=0; i< 100*newNum ; ++i)
+		{
+			realx3 p = pReg.peek();
+			// check if point is inside internal box 
+			if(!internalBox.isInside(p))continue;
+
+			if( !checkForContact(positions, diams, p, d) )
+			{
+				names.push_back(name);
+				positions.push_back(p);
+				diams.push_back(d);
+				numInserted++;
+
+				if( numInserted == newNum ) break;
+				
+				name = mix.getNextShapeName();
+				shapes_.shapeNameToIndex(name, idx);
+				d = shapes_.boundingDiameter(idx);
+			}
+		}
+	}
+	else
+	{
+		real maxDiam = shapes_.maxBoundingSphere();
+		auto minP = pReg.minPoint() - maxDiam;
+		auto maxP = pReg.maxPoint() + maxDiam;
+		auto bDict = dictionary("boxInfo");
+		bDict.add("min", minP);
+		bDict.add("max", maxP);
+		auto selector = pStructSelector::create(
+			"box", 
+			pStruct(), 
+			bDict);
+			
+		auto allPositions = selector().selectedPointPositions();
+		auto allDiameters = selectedFieldVals<real>(selector(), "diameter");
+		auto collCheck = collisionCheck(
+			{minP, maxP},
+			maxDiam,
+			allPositions,
+			allDiameters);
+		
+		uint32 idx;
+		word name = mix.getNextShapeName();
+		shapes_.shapeNameToIndex(name, idx);
+		real d = shapes_.boundingDiameter(idx);
+
+		for(uint32 i=0; i< 100*newNum ; ++i)
+		{
+			realx3 p = pReg.peek();
+			// check if point is inside internal box 
+			if(!internalBox.isInside(p))continue;
+
+			if( collCheck.checkPoint( p, d) )
+			{
+				names.push_back(name);
+				positions.push_back(p);
+				numInserted++;
+				if( numInserted == newNum ) break;
+				
+				// add this new particle to collision check set 
+				allDiameters.push_back(d);
+				allPositions.push_back(p);
+				collCheck.mapLastAddedParticle();
+								
+				// obtain next shape name and diameter 
+				name = mix.getNextShapeName();
+				shapes_.shapeNameToIndex(name, idx);
+				d = shapes_.boundingDiameter(idx);
+			}
+			
+		}
+	}*/

src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp

@@ -21,9 +21,15 @@ Licence:
 #ifndef __InsertionRegion_hpp__
 #define __InsertionRegion_hpp__
 
+#include "insertionRegion.hpp"
 #include "dictionary.hpp"
 #include "pointStructure.hpp"
-#include "insertionRegion.hpp"
+#include "insertion.hpp"
+#include "collisionCheck.hpp"
+#include "pStructSelector.hpp"
+#include "fieldSelector.hpp"
+
+
 
 namespace pFlow
 {

src/Particles/Insertion/collisionCheck/collisionCheck.cpp

@@ -49,16 +49,29 @@ pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
     if(!searchBox_.isInside(p)) return false;
 
     const auto ind = pointIndex(p);
-    
-    uint32 n = head_(ind.x(), ind.y(), ind.z());
-    while( n != -1)
+    const auto startInd = max(ind - 1 ,int32x3(0));
+    const auto endInd = min( ind+1 ,nCells_-1);
+
+    for(int32 i=startInd.x(); i<=endInd.x(); i++)
     {
-        if( (position_[n]-p).length() - 0.5*(diameters_[n]+d )<0.0 )
+        for(int32 j=startInd.y(); j<=endInd.y(); j++)
         {
-            return false;
+            for(int32 k=startInd.z(); k<=endInd.z(); k++)
+            {
+                uint32 n = head_(i, j, k);
+    
+                while( n != -1)
+                {
+                    if( ((position_[n]-p).length() - 0.5*(diameters_[n]+d )) <= 0.0 )
+                    {
+                        return false;
+                    }
+                    n = next_[n];
+                }
+            }
         }
-        n = next_[n];
     }
+    
 	return true;
 }
 

src/Particles/Insertion/insertion/insertion.cpp

@@ -30,13 +30,42 @@ pFlow::insertion::insertion(particles& prtcl)
         insertionFile__,
         "",
         objectFile::READ_IF_PRESENT,
-        objectFile::WRITE_NEVER
+        objectFile::WRITE_ALWAYS
       ),
-      &prtcl.control().caseSetup()
+      &prtcl.time()
     ),
     particles_(prtcl)
 {
-	readInsertionDict(*this);
+	// this means that insertion file exist in time folder 
+	if( IOobject::implyRead() )
+	{
+		readFromFile_ = true;
+		
+	} // look inside the caseSetup folder if it exist
+	else
+	{
+		// read dictionary from caseSetup folder 
+		fileDictionary caseFile(
+			objectFile(
+				insertionFile__,
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_NEVER),
+			&prtcl.control().caseSetup());
+		
+		// check if read happened
+		if(caseFile.implyRead())
+		{
+			readFromFile_ = true;
+			// assign it to this dictionary 
+			fileDictionary::dictionary::operator=(caseFile);
+		}
+	}
+
+	if( readFromFile_)
+	{
+		readInsertionDict();
+	}
 }
 
 const pFlow::pointStructure&
@@ -45,27 +74,18 @@ pFlow::insertion::pStruct() const
 	return particles_.pStruct();
 }
 
-bool
-pFlow::insertion::read(iIstream& is, const IOPattern& iop)
-{
-	if (fileDictionary::read(is, iop))
-	{
-		readFromFile_ = true;
-		return true;
-	}
-	else
-	{
-		return false;
-	}
-}
+
+
 
 bool
-pFlow::insertion::readInsertionDict(const dictionary& dict)
+pFlow::insertion::readInsertionDict()
 {
-	active_ = dict.getVal<Logical>("active");
+	active_ = getVal<Logical>("active");
 
 	if (active_)
 	{
+    	checkForCollision_ = getVal<Logical>("checkForCollision");
+
 		REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
 		          << " in the simulation." << END_REPORT;
 	}
@@ -78,26 +98,56 @@ pFlow::insertion::readInsertionDict(const dictionary& dict)
 	return true;
 }
 
-/*
-bool pFlow::insertion::writeInsertionDict
-(
-                dictionary& dict
-)const
+bool
+pFlow::insertion::writeInsertionDict(dictionary& dict) const
 {
-                if(!dict.add("active", active_) )
-                {
-                                fatalErrorInFunction<<
-                                "  error in writing active to dictionary
-"<<dict.globalName()<<endl; return false;
-                }
+	if (!dict.add("active", active_))
+	{
+		fatalErrorInFunction <<"Error in writing active to dictionary "
+		                        <<dict.globalName()<<endl; 
+		return false;
+	}
 
-                if(!dict.add("checkForCollision", checkForCollision_) )
-                {
-                                fatalErrorInFunction<<
-                                "  error in writing checkForCollision to
-dictionary
-"<<dict.globalName()<<endl; return false;
-                }
+	if (!dict.add("checkForCollision", checkForCollision_))
+	{
+		fatalErrorInFunction << 
+		"Error in writing checkForCollision to dictionary "<<
+		dict.globalName()<<endl; 
+		return false;
+	}
 
-                return true;
+	return true;
+}
+
+
+bool pFlow::insertion::write(iOstream & os, const IOPattern & iop) const
+{
+	dictionary newDict(fileDictionary::dictionary::name(), true);
+	if( iop.thisProcWriteData() )
+	{
+		if( !writeInsertionDict(newDict) || 
+			!newDict.write(os))
+		{
+			fatalErrorInFunction<<
+			" error in writing to dictionary "<< newDict.globalName()<<endl;
+			return false;
+		}
+	}
+
+	return true;
+}
+
+
+/*bool
+pFlow::insertion::read(iIstream& is, const IOPattern& iop)
+{
+	if (fileDictionary::read(is, iop))
+	{
+		readFromFile_ = true;
+		return true;
+	}
+	else
+	{
+		return false;
+	}
 }*/

src/Particles/Insertion/insertion/insertion.hpp

@@ -46,16 +46,22 @@ private:
 	/// to be inserted due to collision
 	Logical    increaseVelocity_ = "No";
 
+	word 		diameterName_ = "diameter";
+
+	word 		velocityName_ = "velocity";
+
 	/// Ref to particles
 	particles& particles_;
 
 	bool       readFromFile_ = false;
 
 	/// Read from dictionary
-	bool       readInsertionDict(const dictionary& dict);
+	bool readInsertionDict();
+
+protected:
 
 	/// Write to dictionary
-	// bool writeInsertionDict(dictionary& dict)const;
+	virtual bool writeInsertionDict(dictionary& dict)const;
 
 public:
 
@@ -71,7 +77,7 @@ public:
 	/// is Insertion active
 	inline bool isActive() const
 	{
-		return active_();
+		return readFromFile_ && active_();
 	}
 
 	inline bool checkForCollision() const
@@ -96,14 +102,24 @@ public:
 		return readFromFile_;
 	}
 
+	const word& diameterName()const
+	{
+		return diameterName_;
+	}
+
+	const word& velocityName()const
+	{
+		return velocityName_;
+	}
+
 	/// read from stream
-	bool read(iIstream& is, const IOPattern& iop) override;
+	//bool read(iIstream& is, const IOPattern& iop) override;
 
 	/*/// read from iIstream
-	virtual bool read(iIstream& is) = 0;
+	virtual bool read(iIstream& is) = 0;*/
 
 	/// write to iOstream
-	virtual bool write(iOstream& os)const = 0;*/
+	bool write(iOstream& os, const IOPattern& iop)const override ;
 };
 
 }

src/Particles/Insertion/insertionRegion/insertionRegion.cpp

@@ -101,9 +101,15 @@ pFlow::insertionRegion::readInsertionRegion(const dictionary& dict)
 bool
 pFlow::insertionRegion::writeInsertionRegion(dictionary& dict) const
 {
-	if (!dict.add("type", type_))
+	
+	if(!dict.add("rate", rate_))
 		return false;
 
+	if(!tControl_.write(dict))
+		return false;
+
+	if (!dict.add("regionType", type_))
+		return false;
 	if (pRegion_)
 	{
 		auto& prDict = dict.subDictOrCreate(type_ + "Info");
@@ -118,11 +124,11 @@ pFlow::insertionRegion::writeInsertionRegion(dictionary& dict) const
 			return false;
 	}
 
-	/*if(setFields_)
+	if(setFieldDict_)
 	{
-	    auto& sfDict = dict.subDictOrCreate("setFields");
-	    setFields_().write(sfDict);
-	}*/
+	    if(!dict.addDict("setFields", setFieldDict_()))
+			return false;	    
+	}
 
 	return true;
 }

src/Particles/Insertion/insertionRegion/insertionRegion.hpp

@@ -159,6 +159,11 @@ public:
 		return type_;
 	}
 
+	const auto& dict()const
+	{
+		return dict_;
+	}
+
 	const auto& Insertion() const
 	{
 		return insertion_;
@@ -206,16 +211,13 @@ public:
 	        return false;
 
 	    return readInsertionRegion(dict);
-	}
+	}*/
 
 	/// write to dictionary
 	bool write(dictionary& dict) const
 	{
-	    if (!timeFlowControl::write(dict))
-	        return false;
-
 	    return writeInsertionRegion(dict);
-	}*/
+	}
 };
 
 } // pFlow

src/Particles/Insertion/shapeMixture/shapeMixture.hpp

@@ -58,7 +58,9 @@ private:
 	uint32Vector   	numberInserted_{"numberInserted"};
 
 	/// Current number of inserted
-	uint32Vector      current_{"currentInserted"};
+	uint32Vector    current_{"currentInserted"};
+
+	uint32 			lastPeaked_ = 0;
 
 public:
 

src/Particles/particles/particleIdHandler/particleIdHandler.cpp

@@ -61,7 +61,7 @@ pFlow::particleIdHandler::hearChanges(
 		auto idRange = getIdRange(numNew);
 		uint32Vector newId("newId",numNew,numNew,RESERVE());
 		fillSequence(newId, idRange.first);
-		output<< "id "<< idRange<<endl;
+		
 		return this->field().insertSetElement(indices, newId);
 	}
 	else

src/phasicFlow/containers/pointField/internalField/internalField.cpp

@@ -26,8 +26,7 @@ bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
 	const auto& indices = varList.getObject<uint32IndexContainer>(
 		eventName);
 	bool success = false;
-	output<<"insert for field "<< name()<<endl;
-
+	
 	if(varList.contains(name()))
 	{
 		// a single value is assigned

src/phasicFlow/dictionary/entry/dataEntry.cpp

@@ -121,6 +121,9 @@ bool pFlow::dataEntry::writeDataEntry
 			}
 			else
 			{
+				
+				if(tok != tokens.cbegin())os<<spaceToken();
+
 				os << *tok;	
 			}
 			lastPuncBegin = true;
@@ -141,6 +144,7 @@ bool pFlow::dataEntry::writeDataEntry
 		}
 		else
 		{
+			
 			if(!lastPuncBegin&&applySpace) os<<spaceToken();
 			os << *tok;
 			lastPuncEnd = false;

src/phasicFlow/repository/Time/baseTimeControl.cpp

@@ -27,6 +27,8 @@ pFlow::baseTimeControl::baseTimeControl
     const word& intervalPrefix,
 	real defStartTime
 )
+:
+	intervalPrefix_(intervalPrefix)
 {
 	auto tControl = dict.getVal<word>("timeControl");
 	if(tControl == "timeStep")
@@ -162,3 +164,32 @@ pFlow::baseTimeControl::iInterval() const
 	fatalExit;
 	return 0;
 }
+
+bool
+pFlow::baseTimeControl::write(dictionary& dict) const
+{
+	if(isTimeStep_)
+	{
+		dict.add("timeControl", "timeStep");
+	}
+	else
+	{
+		dict.add("timeControl", "runTime");
+	}
+
+	word intervalWord = intervalPrefix_.size()==0? word("interval"): intervalPrefix_+"Interval";
+	
+	if(!isTimeStep_)
+	{
+		dict.add(intervalWord,rRange_.stride());
+		dict.add("startTime",rRange_.begin());
+		dict.add("endTime", rRange_.end());
+	}
+	else
+	{
+		dict.add(intervalWord,iRange_.stride());
+		dict.add("startTime",iRange_.begin());
+		dict.add("endTime", iRange_.end());
+	}
+	return true;
+}

src/phasicFlow/repository/Time/baseTimeControl.hpp

@@ -1,20 +1,20 @@
 /*------------------------------- phasicFlow ---------------------------------
-			O        C enter of
-		 O O       E ngineering and
-		O   O      M ultiscale modeling of
-	 OOOOOOO     F luid flow       
+            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
 ------------------------------------------------------------------------------
-	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.
+    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 __baseTimeControl_hpp__
@@ -23,50 +23,51 @@ Licence:
 #include "dictionary.hpp"
 #include "ranges.hpp"
 
-
 namespace pFlow
 {
 
 class baseTimeControl
 {
 private:
-	 
-	bool 	isTimeStep_;
+
+	bool              isTimeStep_;
 
 	int32StridedRagne iRange_;
 
 	realStridedRange  rRange_;
-	
+
+	const word        intervalPrefix_;
 
 public:
 
 	baseTimeControl(
-		const dictionary& dict, 
-		const word& intervalPrefix = "",
-		real defStartTime = 0.0);
-	
-	inline
-	bool isTimeStep()const
+	  const dictionary& dict,
+	  const word&       intervalPrefix = "",
+	  real              defStartTime   = 0.0
+	);
+
+	inline bool isTimeStep() const
 	{
 		return isTimeStep_;
 	}
 
-	bool timeEvent(uint32 iter, real t, real dt)const;
+	bool  timeEvent(uint32 iter, real t, real dt) const;
 
-	bool isInRange(uint32 iter, real t, real dt)const;
+	bool  isInRange(uint32 iter, real t, real dt) const;
 
-	real startTime()const;
+	real  startTime() const;
 
-	real endTime()const;
+	real  endTime() const;
 
-	real rInterval()const;
+	real  rInterval() const;
 
-	int32 startIter()const;
+	int32 startIter() const;
 
-	int32 endIter()const;
+	int32 endIter() const;
 
-	int32 iInterval()const;
+	int32 iInterval() const;
 
+	bool  write(dictionary& dict) const;
 };
 
 }

src/phasicFlow/structuredData/cylinder/cylinder.cpp

@@ -21,11 +21,13 @@ Licence:
 
 #include "cylinder.hpp"
 #include "zAxis.hpp"
+#include "streams.hpp"
 
 FUNCTION_H
 bool pFlow::cylinder::calculateParams()
 {
 
+	WARNING<<"Use of cylinder requires modifications to zAxis"<<END_WARNING;
 	auto p1p2 = p2_ - p1_;
 
 	if( p1p2.length() > smallValue )

src/phasicFlow/structuredData/pointStructure/selectors/pStructSelector/fieldSelector.hpp

@@ -33,7 +33,7 @@ Vector<T> selectedFieldVals(const pStructSelector& selector, const word& name)
     {
         fatalErrorInFunction<<
 		"Type of field "<< name << " in time repository does not match with"<<
-        FieldType::TYEPANME()<<endl;
+        FieldType::TYPENAME()<<endl;
 		fatalExit;
     }
 

src/phasicFlow/structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp

@@ -32,6 +32,16 @@ pFlow::pStructSelector::pStructSelector
 	pStruct_(pStruct)
 {}
 
+pFlow::pStructSelector::pStructSelector(
+  const word&           ,
+  const pointStructure& pStruct,
+  const dictionary&     
+)
+:
+	pStruct_(pStruct)
+{
+}
+
 const pFlow::pointStructure& pFlow::pStructSelector::pStruct()const
 {
 	return pStruct_;
@@ -63,7 +73,6 @@ pFlow::pStructSelector::create(
   const dictionary&     dict
 )
 {
-
 	word selectorMethod = angleBracketsNames("selector", dict.getVal<word>("selector"));
 
 	if( dictionaryvCtorSelector_.search(selectorMethod) )
@@ -74,7 +83,7 @@ pFlow::pStructSelector::create(
 	{
 		printKeys
 		( 
-			fatalError << "Ctor Selector "<< selectorMethod << " dose not exist. \n"
+			fatalError << "Ctor Selector "<< selectorMethod << " does not exist. \n"
 			<<"Avaiable ones are: \n\n"
 			,
 			dictionaryvCtorSelector_
@@ -83,3 +92,31 @@ pFlow::pStructSelector::create(
 	}
 	return nullptr;
 }
+
+pFlow::uniquePtr<pFlow::pStructSelector>
+pFlow::pStructSelector::create(
+  const word&           type,
+  const pointStructure& pStruct,
+  const dictionary&     dict
+)
+{
+	word selectorMethod = angleBracketsNames("selector", type);
+
+	if( wordvCtorSelector_.search(selectorMethod) )
+	{
+		return wordvCtorSelector_[selectorMethod] (type, pStruct, dict);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< selectorMethod << " does not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		fatalExit;
+	}
+	return nullptr;
+
+}

src/phasicFlow/structuredData/pointStructure/selectors/pStructSelector/pStructSelector.hpp

@@ -44,9 +44,16 @@ public:
 	// - type info
 	TypeInfo("pStructSelector");
 
-
+		/// the dictionary contains the selector keyword and another dictionary which is 
+		/// used for creating selector 
 		pStructSelector(const pointStructure& pStruct, const dictionary& UNUSED(dict));
 
+		/// construct using selector type and a dictionary that contains info of selector 
+		pStructSelector(
+			const word& type, 
+			const pointStructure& pStruct, 
+			const dictionary& UNUSED(dict));
+
 		create_vCtor
 		(
 			pStructSelector,
@@ -55,6 +62,18 @@ public:
 			(pStruct, dict)
 		);
 
+		create_vCtor
+		(
+			pStructSelector,
+			word,
+			(
+				const word& type, 
+				const pointStructure& pStruct, 
+				const dictionary& dict
+			),
+			(type, pStruct, dict)
+		);
+
 		virtual ~pStructSelector() = default;
 
 	//// - Methods
@@ -71,6 +90,11 @@ public:
 	static
 	uniquePtr<pStructSelector> create(const pointStructure& pStruct, const dictionary& dict);
 
+	static
+	uniquePtr<pStructSelector> create(
+		const word& type, 
+		const pointStructure& pStruct, 
+		const dictionary& dict);
 };
 
 

src/phasicFlow/structuredData/pointStructure/selectors/selectorGeometric/selectorGeometric.cpp

@@ -1,20 +1,18 @@
 
 
-
 template<typename GeomType>
 bool
 pFlow::selectorGeometric<GeomType>::selectPoints()
 {
-	
 	selectedPoints_.clear();
 
-	auto maskH = pStruct().activePointsMaskHost();
+	auto       maskH  = pStruct().activePointsMaskHost();
 	const auto aRange = maskH.activeRange();
-	auto pPos = pStruct().pointPositionHost();
-	
-	for(uint32 i=aRange.start(); i<aRange.end();i++)
+	auto       pPos   = pStruct().pointPositionHost();
+
+	for (uint32 i = aRange.start(); i < aRange.end(); i++)
 	{
-		if(maskH.isActive(i)&& pRegion_.isInside( pPos[i] ))
+		if (maskH.isActive(i) && pRegion_.isInside(pPos[i]))
 		{
 			selectedPoints_.push_back(i);
 		}
@@ -32,10 +30,24 @@ pFlow::selectorGeometric<GeomType>::selectorGeometric(
     type_(dict.getVal<word>("selector")),
     pRegion_(type_, dict.subDict(type_ + "Info"))
 {
-
-	if(!selectPoints())
+	if (!selectPoints())
 	{
 		fatalExit;
 	}
 }
 
+template<typename GeomType>
+pFlow::selectorGeometric<GeomType>::selectorGeometric(
+  const word&           type,
+  const pointStructure& pStruct,
+  const dictionary&     dict
+)
+  : pStructSelector(type, pStruct, dict),
+    type_(type),
+    pRegion_(type_, dict)
+{
+	if (!selectPoints())
+	{
+		fatalExit;
+	}
+}

src/phasicFlow/structuredData/pointStructure/selectors/selectorGeometric/selectorGeometric.hpp

@@ -50,6 +50,11 @@ public:
 
 	selectorGeometric(const pointStructure& pStruct, const dictionary& dict);
 
+	selectorGeometric(
+		const word& type, 
+		const pointStructure& pStruct, 
+		const dictionary& dict);
+		
 	add_vCtor
     (
         pStructSelector, 
@@ -57,6 +62,13 @@ public:
         dictionary
     );
 
+	add_vCtor
+    (
+        pStructSelector, 
+        selectorGeometric, 
+        word
+    );
+
 	~selectorGeometric() override = default;
 
 	const uint32Vector& selectedPoints()const override

src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.cpp

@@ -2,56 +2,54 @@
       O        C enter of
      O O       E ngineering and
     O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
+   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 
+  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 
+  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 "selectorRandomPoints.hpp"
+#include "Set.hpp"
 #include "dictionary.hpp"
 #include "uniformRandomUint32.hpp"
-#include "Set.hpp"
 
-
-bool pFlow::selectorRandomPoints::selectAllPointsInRange()
+bool
+pFlow::selectorRandomPoints::selectAllPointsInRange()
 {
 	// to reduct allocations
-	uint32 maxNum = number_+1;
-	
-	selectedPoints_.reserve	(maxNum);
+	uint32 maxNum = number_ + 1;
+
+	selectedPoints_.reserve(maxNum);
 
 	selectedPoints_.clear();
-	
-	uniformRandomUint32 intRand (begin_, end_);
 
+	uniformRandomUint32 intRand(begin_, end_);
 
-	uint32 n = 0;
-	uint32 iter = 0;
-	bool finished = false;
+	uint32              n        = 0;
+	uint32              iter     = 0;
+	bool                finished = false;
 
-	Set<uint32> selctedIndices;
+	Set<uint32>         selctedIndices;
 
-	while( iter < number_*100)
+	while (iter < number_ * 100)
 	{
 		uint32 newInd = intRand.randomNumber();
 
-		if( auto [it, inserted] = selctedIndices.insert(newInd); inserted )
+		if (auto [it, inserted] = selctedIndices.insert(newInd); inserted)
 		{
 			n++;
 
-			if(n == number_)
+			if (n == number_)
 			{
 				finished = true;
 				break;
@@ -60,62 +58,61 @@ bool pFlow::selectorRandomPoints::selectAllPointsInRange()
 		iter++;
 	}
 
-	
-	if(finished)
+	if (finished)
 	{
-		for(auto& ind:selctedIndices)
+		for (auto& ind : selctedIndices)
 		{
 			selectedPoints_.push_back(ind);
 		}
 
 		return true;
-
-	}else
-	{
-		fatalErrorInFunction<< "Could not find random indices in the range."<<endl;
-		return false;	
 	}
-
-	
+	else
+	{
+		fatalErrorInFunction << "Could not find random indices in the range."
+		                     << endl;
+		return false;
+	}
 }
-	
-pFlow::selectorRandomPoints::selectorRandomPoints
-(
-	const pointStructure& pStruct,
-	const dictionary& dict
+
+pFlow::selectorRandomPoints::selectorRandomPoints(
+  const pointStructure& pStruct,
+  const dictionary&     dict
 )
-:
-	pStructSelector
-	(
-		pStruct, dict
-	),
-	begin_
-	(
-		dict.subDict("randomPointsInfo").getVal<uint32>("begin")
-	),
-	end_
-	(
-		dict.subDict("randomPointsInfo").getValOrSet("end", pStruct.size())
-	),
-	number_
-	(
-		dict.subDict("randomPointsInfo").getValOrSet("number", 1)
-	)
+  : selectorRandomPoints(
+      "randomPoints",
+      pStruct,
+      dict.subDict("randomPointsInfo")
+    )
 {
-	begin_ 	= max(begin_,1u);
-	end_ 	= min(end_, static_cast<uint32>(pStruct.size()));
-	number_ = max(number_,0u);
-	if(end_-begin_ < number_)
+}
+
+pFlow::selectorRandomPoints::selectorRandomPoints(
+  const word&           type,
+  const pointStructure& pStruct,
+  const dictionary&     dict
+)
+  : pStructSelector(type, pStruct, dict),
+    begin_(dict.getVal<uint32>("begin")),
+    end_(dict.getValOrSet("end", pStruct.size())),
+    number_(dict.getValOrSet("number", 1))
+{
+	begin_  = max(begin_, 1u);
+	end_    = min(end_, static_cast<uint32>(pStruct.size()));
+	number_ = max(number_, 0u);
+
+	if (end_ - begin_ < number_)
 	{
+		warningInFunction
+		  << "In dictionary " << dict.globalName()
+		  << " number is greater than the interval defined by begine and end ["
+		  << begin_ << " " << end_ << "), resetting it to " << end_ - begin_
+		  << endl;
 
-		warningInFunction<< "In dictionary " << dict.globalName()<<
-		" number is greater than the interval defined by begine and end ["<<
-		begin_<<" "<<end_<<"), resetting it to "<<end_-begin_<<endl;
-
-		number_ = end_-begin_;
+		number_ = end_ - begin_;
 	}
 
-	if(!selectAllPointsInRange())
+	if (!selectAllPointsInRange())
 	{
 		fatalExit;
 	}

src/phasicFlow/structuredData/pointStructure/selectors/selectorRandomPoints/selectorRandomPoints.hpp

@@ -58,6 +58,11 @@ public:
 
 		selectorRandomPoints(const pointStructure& pStruct, const dictionary& dict);
 
+		selectorRandomPoints(
+			const word& type, 
+			const pointStructure& pStruct, 
+			const dictionary& UNUSED(dict));
+
 		add_vCtor
 		(
 			pStructSelector,
@@ -65,6 +70,12 @@ public:
 			dictionary
 		);
 
+		add_vCtor
+		(
+			pStructSelector,
+			selectorRandomPoints,
+			word
+		);
 		~selectorRandomPoints() final = default;
 
 	//// - Methods

src/phasicFlow/structuredData/pointStructure/selectors/selectorStridedRange/selectorStridedRange.cpp

@@ -2,67 +2,66 @@
       O        C enter of
      O O       E ngineering and
     O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
+   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 
+  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 
+  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 "selectorStridedRange.hpp"
 #include "dictionary.hpp"
 
-void pFlow::selectorStridedRange::selectAllPointsInRange()
+void
+pFlow::selectorStridedRange::selectAllPointsInRange()
 {
 	// to reduct allocations
-	uint32 maxNum = (end_ - begin_)/stride_+2;
-	
-	selectedPoints_.reserve	(maxNum);
+	uint32 maxNum = (end_ - begin_) / stride_ + 2;
+
+	selectedPoints_.reserve(maxNum);
 
 	selectedPoints_.clear();
-		
-	for(uint32 i = begin_; i<end_; i+= stride_)
+
+	for (uint32 i = begin_; i < end_; i += stride_)
 	{
 		selectedPoints_.push_back(i);
 	}
 }
 
-pFlow::selectorStridedRange::selectorStridedRange
-(
-	const pointStructure& pStruct,
-	const dictionary& dict
+pFlow::selectorStridedRange::selectorStridedRange(
+  const pointStructure& pStruct,
+  const dictionary&     dict
 )
-:
-	pStructSelector
-	(
-		pStruct, dict
-	),
-	begin_
-	(
-		dict.subDict("stridedRangeInfo").getVal<uint32>("begin")
-	),
-	end_
-	(
-		dict.subDict("stridedRangeInfo").getValOrSet<uint32>("end", pStruct.size())
-	),
-	stride_
-	(
-		dict.subDict("stridedRangeInfo").getValOrSet<uint32>("stride", 1u)
-	)
+  : selectorStridedRange(
+      "stridedRange",
+      pStruct,
+      dict.subDict("stridedRangeInfo")
+    )
 {
-	begin_ 	= max(begin_,1u);
-	end_ 	= min(end_, static_cast<uint32>(pStruct.size()));
-	stride_ = max(stride_,1u);
+}
+
+pFlow::selectorStridedRange::selectorStridedRange(
+  const word&           type,
+  const pointStructure& pStruct,
+  const dictionary&     dict
+)
+  : pStructSelector(type, pStruct, dict),
+    begin_(dict.getVal<uint32>("begin")),
+    end_(dict.getValOrSet<uint32>("end", pStruct.size())),
+    stride_(dict.getValOrSet<uint32>("stride", 1u))
+{
+	begin_  = max(begin_, 1u);
+	end_    = min(end_, static_cast<uint32>(pStruct.size()));
+	stride_ = max(stride_, 1u);
 
 	selectAllPointsInRange();
 }

src/phasicFlow/structuredData/pointStructure/selectors/selectorStridedRange/selectorStridedRange.hpp

@@ -58,6 +58,11 @@ public:
 
 		selectorStridedRange(const pointStructure& pStruct, const dictionary& dict);
 
+		selectorStridedRange(
+			const word& type, 
+			const pointStructure& pStruct, 
+			const dictionary& dict);
+		
 		add_vCtor
 		(
 			pStructSelector,
@@ -65,6 +70,13 @@ public:
 			dictionary
 		);
 
+		add_vCtor
+		(
+			pStructSelector,
+			selectorStridedRange,
+			word
+		);
+
 		~selectorStridedRange() final = default;
 
 	//// - Methods

src/phasicFlow/types/triple/triple.hpp

@@ -50,17 +50,17 @@ struct triple
 	T y_;
 	T z_;
 
+	/// Type info for triple
 	TripleTypeInfoNV(T);
+	
 	//// Constructors
 
 	/// Initilize to zero
 	INLINE_FUNCTION_HD
-	triple()
-	  : x_(),
-	    y_(),
-	    z_()
-	{
-	}
+	triple() = default;
+	
+	INLINE_FUNCTION_HD
+	~triple() = default;
 
 	/// Construct from x, y, z
 	INLINE_FUNCTION_HD