Particle insertion is added with anyList

- collision check is not active yet.
- variable velocity is not active yet.
- events and messages are not active yet.

solvers/sphereGranFlow/createDEMComponents.hpp

@@ -36,10 +36,9 @@ REPORT(0)<<"\nCreating particle insertion object . . ."<<END_REPORT;
 		sphParticles.shapes()
 	);*/
 
-/*auto sphInsertion =  sphereInsertion(
+auto sphInsertion =  sphereInsertion(
-	Control.caseSetup().path()+insertionFile__, 
 	sphParticles, 
-	sphParticles.shapes());*/
+	sphParticles.spheres());
 
 REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<END_REPORT;
 auto interactionPtr = interaction::create(

solvers/sphereGranFlow/sphereGranFlow.cpp

@@ -38,7 +38,7 @@ Licence:
 #include "geometry.hpp"
 #include "sphereParticles.hpp"
 #include "interaction.hpp"
-//#include "Insertions.hpp"
+#include "Insertions.hpp"
 
 
 
@@ -80,14 +80,15 @@ initialize_pFlowProcessors();
 	do 
 	{
 		
-		/*if(! sphInsertion.insertParticles( 
+		if(! sphInsertion.insertParticles( 
+			Control.time().currentIter(),
 			Control.time().currentTime(),
 			Control.time().dt()	) )
 		{
 			fatalError<<
 			"particle insertion failed in sphereDFlow solver.\n";
 			return 1;
-		}*/
+		}
 		 
 		// set force to zero
 		surfGeometry.beforeIteration();

src/Integration/AdamsBashforth2/AdamsBashforth2.hpp

@@ -48,8 +48,8 @@ private:
 	
 public:
 
-	/// Type info
+	/// Class info
-	TypeInfo("AdamsBashforth2");
+	ClassInfo("AdamsBashforth2");
 
 	// - Constructors
 
@@ -95,8 +95,16 @@ public:
 		bool correct(
 			real dt, 
 			realx3Field_D& y, 
-			realx3PointField_D& dy);
+			realx3PointField_D& dy) final;
 		
+		/*bool hearChanges
+		(
+			real t,
+			real dt,
+			uint32 iter,
+			const message& msg, 
+			const anyList& varList
+		) override;*/
 
 		bool setInitialVals(
 			const int32IndexContainer& newIndices,

src/Particles/CMakeLists.txt

@@ -1,19 +1,18 @@
 
 set(SourceFiles 
 dynamicPointStructure/dynamicPointStructure.cpp
-particles/baseShapeNames.cpp
+particles/shape/baseShapeNames.cpp
-particles/shape.cpp
+particles/shape/shape.cpp
 particles/particles.cpp
-particles/particleIdHandler.cpp
+particles/particleIdHandler/particleIdHandler.cpp
-particles/regularParticleHandler.cpp
+particles/regularParticleIdHandler/regularParticleIdHandler.cpp
 SphereParticles/sphereShape/sphereShape.cpp
 SphereParticles/sphereParticles/sphereParticles.cpp
 SphereParticles/sphereParticles/sphereParticlesKernels.cpp
-#Insertion/shapeMixture/shapeMixture.cpp
+Insertion/insertionRegion/insertionRegion.cpp
-#Insertion/insertion/insertion.cpp
+Insertion/insertion/insertion.cpp
-#Insertion/Insertion/Insertions.cpp
+Insertion/shapeMixture/shapeMixture.cpp
-#Insertion/insertionRegion/insertionRegion.cpp
+Insertion/Insertion/Insertions.cpp
-#Insertion/insertionRegion/timeFlowControl.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Integration Property)

src/Particles/Insertion/Insertion/Insertion.cpp

@@ -18,7 +18,7 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-template<typename ShapeType>
+/*template<typename ShapeType>
 bool pFlow::Insertion<ShapeType>::readInsertionDict
 (
 	const dictionary& dict
@@ -65,6 +65,33 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
 	}
 
 	return true;
+}*/
+
+template<typename ShapeType>
+bool
+pFlow::Insertion<ShapeType>::setInsertionRegions()
+{
+	regions_.clear();
+	
+	if( !this->isActive() )
+	{
+		return true;
+	}
+
+	wordList regionDicNames = this->dictionaryKeywords();
+
+	for(const auto& name:regionDicNames)
+	{
+		REPORT(2)<<"reading insertion region "<< Green_Text(name)<<END_REPORT;
+		
+		regions_.push_back(makeUnique<InsertionRegion<ShapeType>>(
+			name,
+			*this, 
+			shapes_));
+	}
+
+	return true;
+
 }
 
 template<typename ShapeType>
@@ -75,11 +102,14 @@ pFlow::Insertion<ShapeType>::Insertion(
 	insertion(prtcl),
 	shapes_(shapes)
 {
-
+	if(!setInsertionRegions())
-	
+	{
+		fatalErrorInFunction;
+		fatalExit;
+	}	
 }
 
-template<typename ShapeType>
+/*template<typename ShapeType>
 pFlow::Insertion<ShapeType>::Insertion(
 	fileSystem file, 
 	particles& prtcl, 
@@ -95,13 +125,14 @@ pFlow::Insertion<ShapeType>::Insertion(
 		file<<endl;
 		fatalExit;
 	}
-}
+}*/
 
 
 template<typename ShapeType>
 bool pFlow::Insertion<ShapeType>::insertParticles
 (
-	real currentTime,
+	uint32 iter, 
+	real t,
 	real dt
 )
 {
@@ -112,21 +143,21 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 	{
 		bool insertionOccured = false;
 		auto& rgn = regions_[i];
-		if( rgn.insertionTime(currentTime, dt) )
+		if( rgn.insertionTime(iter, t, dt) )
 		{
 			
 			realx3Vector pos;
 			wordVector shapes;
-			if( rgn.insertParticles(currentTime, dt, shapes, pos, insertionOccured) )
+			if( rgn.insertParticles(iter, t, dt, shapes, pos, insertionOccured) )
 			{
 
 				if(insertionOccured)
 				{
-					REPORT(0)<<"\nParticle insertion from "<< greenText(rgn.name())<<endREPORT;
+					REPORT(0)<<"\nParticle insertion from "<< Green_Text(rgn.name())<<END_REPORT;
-					REPORT(1)<< cyanText(pos.size()) << " new particles is being inserted at Time: "<<
+					REPORT(1)<< Cyan_Text(pos.size()) << " new particles is being inserted at Time: "<<
-						cyanText(currentTime) <<" s."<<endREPORT;
+						Cyan_Text(t) <<" s."<<END_REPORT;
 						
-					if(!particles_.insertParticles(pos, shapes, rgn.setFields()))
+					if(!Particles().insertParticles(pos, shapes, rgn.setFieldList()))
 					{
 						fatalErrorInFunction<<
 						" Cannot add "<< pos.size() << " particles from region "<< rgn.name() << 
@@ -134,7 +165,7 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 						return false;
 					}
 					REPORT(1)<<"Total number of particles inserted from this region is "<<
-						cyanText(rgn.totalInserted())<<'\n'<<endREPORT;
+						Cyan_Text(rgn.totalInserted())<<'\n'<<END_REPORT;
 				}
 				else
 				{
@@ -146,9 +177,15 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 			{
 				if(insertionOccured)
 				{
-					yWARNING<< "\n fewer number of particles are inserted from region "<< rgn.name() <<
+					WARNING<< "\n fewer number of particles are inserted from region "<< rgn.name() <<
-					" than expected. You may stop the simulation to change settings."<<endyWARNING;
+					" than expected. You may stop the simulation to change settings."<<END_WARNING;
-					continue;
+					if(!Particles().insertParticles(pos, shapes, rgn.setFieldList()))
+					{
+						fatalErrorInFunction<<
+						" Cannot add "<< pos.size() << " particles from region "<< rgn.name() << 
+						" to particles. \n";
+						return false;
+					}
 				}
 				else
 				{
@@ -166,7 +203,7 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 }
 
 
-template<typename ShapeType>
+/*template<typename ShapeType>
 bool pFlow::Insertion<ShapeType>::read
 (
 	iIstream& is
@@ -217,4 +254,4 @@ bool pFlow::Insertion<ShapeType>::write
 	}	
 
 	return true;
-}
+}*/

src/Particles/Insertion/Insertion/Insertion.hpp

@@ -18,14 +18,12 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #ifndef __Insertion_hpp__
 #define __Insertion_hpp__
 
-
-#include "insertion.hpp"
-#include "ListPtr.hpp"
 #include "InsertionRegion.hpp"
+#include "ListPtr.hpp"
+#include "insertion.hpp"
 #include "particles.hpp"
 
 namespace pFlow
@@ -37,7 +35,8 @@ namespace pFlow
  *
  * Any number of insertion regions can be defined in a simulation. The
  * data for particle insertion is provided in particleInsertion file, which
- * looks like this. A list of insertion regions (class insertionRegion) can be defined in this file. 
+ * looks like this. A list of insertion regions (class insertionRegion) can be
+defined in this file.
  * For more information see file insertionRegion.hpp.
  * \verbatim
 active  yes;
@@ -54,11 +53,9 @@ region2
  \endverbatim
  */
 template<typename ShapeType>
-class Insertion
+class Insertion : public insertion
-:
-	public insertion
 {
-protected:
+private:
 
 	const ShapeType&                    shapes_;
 
@@ -66,25 +63,25 @@ protected:
 	ListPtr<InsertionRegion<ShapeType>> regions_;
 
 
-	bool readInsertionDict(const dictionary& dict);
+	bool setInsertionRegions();
 
-	bool writeInsertionDict(dictionary& dict)const;
+	// bool readInsertionDict(const dictionary& dict);
+
+	// bool writeInsertionDict(dictionary& dict)const;
 
 public:
 
-	TypeInfoTemplateNV("Insertion",ShapeType);
+	TypeInfoTemplateNV11("Insertion", ShapeType);
 
 	Insertion(particles& prtcl, const ShapeType& shapes);
 
-	Insertion(fileSystem file, particles& prtcl, const ShapeType& shapes);
+	// Insertion(fileSystem file, particles& prtcl, const ShapeType& shapes);
 
+	bool insertParticles(uint32 iter, real t, real dt);
 
-	bool insertParticles(real currentTime, real dt);
+	/*virtual bool read(iIstream& is) override;
-
-	virtual bool read(iIstream& is) override;
-
-	virtual bool write(iOstream& os)const override;
 
+	virtual bool write(iOstream& os)const override;*/
 };
 
 }

src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp

@@ -39,21 +39,21 @@ bool pFlow::InsertionRegion<ShapeType>::checkForContact
 template<typename ShapeType>
 pFlow::InsertionRegion<ShapeType>::InsertionRegion
 (
-	const dictionary& dict,
+	const word& name,
+	const insertion& instn, 
 	const ShapeType& shapes
 )
 :
-	insertionRegion(dict),
+	insertionRegion(name, instn),
 	shapes_(shapes)
-{
+{}
-
-}
 
 
 template<typename ShapeType>
 bool pFlow::InsertionRegion<ShapeType>::insertParticles
 (
-	real currentTime,
+	uint32 iter, 
+	real t,
 	real dt,
 	wordVector& names,
 	realx3Vector& pos,
@@ -62,9 +62,9 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
 {
 	insertionOccured = false;
 
-	if(!insertionTime( currentTime, dt)) return true;
+	if(!insertionTime(iter, t, dt)) return true;
 
-	size_t newNum = numberToBeInserted(currentTime);
+	uint32 newNum = numberToBeInserted(iter, t, dt);
 	
 	if(newNum == 0) return true;
 
@@ -73,43 +73,49 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
 	names.clear();
 	pos.clear();
 
-	realVector diams(newNum, RESERVE());
+	realVector diams("diams", newNum, 0, RESERVE());
 
-	mixture_->getNextShapeNameN(newNum, names);
+	for(uint32 i=0; i<newNum; i++)
-
-	if(!shapes_.shapeToDiameter(names,diams))
 	{
-		fatalErrorInFunction<<
+		uint32 idx;
-		"  error occured in insertion region "<< name() << 
+		shapes_.shapeNameToIndex(names[i], idx);
-		" while converting shapes to diameter. \n";
+		diams[i] = shapes_.boundingDiameter(idx);
-		return false;
 	}
 
 	size_t n = 0;
+	uint32 idx;
+	auto& mix = mixture();
+	auto& pReg = pRegion();
+	word name = mix.getNextShapeName();
+	shapes_.shapeNameToIndex(name, idx);
+	real d = shapes_.boundingDiameter(idx);
 
-	for(label iter=0; iter< 10*newNum ; ++iter)
+	for(uint32 i=0; i< 100*newNum ; ++i)
 	{
-		if( !(n < newNum) ) 
+		
+		realx3 p = pReg.peek();
+		if( !checkForContact(pos, diams, p, d) )
+		{
+			names.push_back(name);
+			pos.push_back(p);
+			diams.push_back(d);
+			n++;
+
+			if( n == newNum ) 
 			{
 				addToNumInserted(newNum);
 				insertionOccured = true;
 				return true;
 			}
-		realx3 p = pRegion_().peek();
+
-		real d = diams[pos.size()];
+			name = mix.getNextShapeName();
-		if( !checkForContact(pos, diams, p, d) )
+			shapes_.shapeNameToIndex(name, idx);
-		{
+			d = shapes_.boundingDiameter(idx);
-			pos.push_back(p);
-			n++;
 		}
 
 	}
 
-	fatalErrorInFunction<<
-	"  Cannot insert "<< newNum << " new particles from region "<< name()<<". \n"
-	"  pFlow could position only "<< n<< " particles in this region. \n";
 	addToNumInserted(n);
-	insertionOccured = false;
+	insertionOccured = true;
 	return false;
-
 }

src/Particles/Insertion/InsertionRegion/InsertionRegion.hpp

@@ -21,9 +21,8 @@ Licence:
 #ifndef __InsertionRegion_hpp__
 #define __InsertionRegion_hpp__
 
-
-#include "insertionRegion.hpp"
 #include "dictionary.hpp"
+#include "insertionRegion.hpp"
 
 namespace pFlow
 {
@@ -34,10 +33,9 @@ namespace pFlow
  */
 template<typename ShapeType>
 class InsertionRegion 
-:
+  : public insertionRegion
-	public insertionRegion
 {
-protected:
+private:
 
 	/// Ref to Shapes
 	const ShapeType& shapes_;
@@ -46,51 +44,34 @@ protected:
 	       const realx3Vector& pos,
 	       const realVector&   diams,
 	       const realx3&       p,
-				const real& d);
+	       const real&         d
+	     );
 
 public:
 
 	/// Type info
-	TypeInfoTemplateNV("insertionRegion", ShapeType);
+	TypeInfoTemplateNV11("insertionRegion", ShapeType);
 
 	// - Constructors
 
 	/// Construct from dictionary
-		InsertionRegion(const dictionary& dict, const ShapeType& shapes);
+	InsertionRegion(
+	  const word&      name,
+	  const insertion& instn,
+	  const ShapeType& shapes
+	);
 
-		/// Copy
+	~InsertionRegion() = default;
-		InsertionRegion(const InsertionRegion<ShapeType>& ) = default;
-
-		/// Move
-		InsertionRegion(InsertionRegion<ShapeType>&&) = default;
-
-		/// Copy assignment
-		InsertionRegion<ShapeType>& operator=(const InsertionRegion<ShapeType>& ) = default;
-
-		/// Copy assignment
-		InsertionRegion<ShapeType>& operator=(InsertionRegion<ShapeType>&&) = default;
-
-		/// Clone
-		auto clone()const
-		{
-			return makeUnique<InsertionRegion<ShapeType>>(*this);
-		}
-
-		/// Clone ptr
-		auto clonePtr()const
-		{
-			return new InsertionRegion<ShapeType>(*this);
-		}
 
 	// - Methods
 
-		/// Insert particles at currentTime
+	/// Insert particles at current time t
 	/// Check if currentTime is the right moment for
 	/// particle insertion. Fill the vectors name, pos and signal
 	/// if particle insertion occured or not.
-		bool insertParticles
+	bool insertParticles(
-		(
+	  uint32        iter,
-			real currentTime,
+	  real          t,
 	  real          dt,
 	  wordVector&   names,
 	  realx3Vector& pos,
@@ -100,14 +81,10 @@ public:
 	// bool read(const dictionary& dict);
 
 	// bool write(dictionary& dict)const;
-
 };
 
-
-
 } // pFlow
 
-
 #include "InsertionRegion.cpp"
 
 #endif

src/Particles/Insertion/insertion/insertion.cpp

@@ -18,31 +18,61 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
-#include "particles.hpp"
-#include "dictionary.hpp"
 #include "insertion.hpp"
+#include "particles.hpp"
 #include "streams.hpp"
+#include "systemControl.hpp"
+#include "vocabs.hpp"
 
-bool pFlow::insertion::readInsertionDict
+pFlow::insertion::insertion(particles& prtcl)
-(
+  : fileDictionary(
-	const dictionary& dict
+      objectFile(
-)
+        insertionFile__,
+        "",
+        objectFile::READ_IF_PRESENT,
+        objectFile::WRITE_NEVER
+      ),
+      &prtcl.control().caseSetup()
+    ),
+    particles_(prtcl)
 {
+	readInsertionDict(*this);
+}
 
+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)
+{
 	active_ = dict.getVal<Logical>("active");
 
 	if (active_)
-		REPORT(1)<< "Particle insertion mechanism is "<<
+	{
-		yellowText("active")<<" in the simulation."<<endREPORT;
+		REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
+		          << " in the simulation." << END_REPORT;
+	}
 	else
-		REPORT(1)<< "Particle insertion mechanism is "<<
+	{
-		yellowText("not active")<<" in the simulation."<<endREPORT;
+		REPORT(1) << "Particle insertion mechanism is "
-
+		          << Yellow_Text("not active") << " in the simulation."
-	
+		          << END_REPORT;
+	}
 	return true;
 }
 
+/*
 bool pFlow::insertion::writeInsertionDict
 (
                 dictionary& dict
@@ -51,24 +81,17 @@ bool pFlow::insertion::writeInsertionDict
                 if(!dict.add("active", active_) )
                 {
                                 fatalErrorInFunction<<
-		"  error in writing active to dictionary "<<dict.globalName()<<endl;
+                                "  error in writing active to dictionary
-		return false;
+"<<dict.globalName()<<endl; return false;
                 }
 
                 if(!dict.add("checkForCollision", checkForCollision_) )
                 {
                                 fatalErrorInFunction<<
-		"  error in writing checkForCollision to dictionary "<<dict.globalName()<<endl;
+                                "  error in writing checkForCollision to
-		return false;
+dictionary
+"<<dict.globalName()<<endl; return false;
                 }
 
                 return true;
-}
+}*/
-
-pFlow::insertion::insertion
-(
-	particles& prtcl
-)
-:
-	particles_(prtcl)
-{}

src/Particles/Insertion/insertion/insertion.hpp

@@ -17,26 +17,23 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-
 #ifndef __insertion_hpp__
 #define __insertion_hpp__
 
-#include "types.hpp"
+#include "fileDictionary.hpp"
-#include "virtualConstructor.hpp"
 
 namespace pFlow
 {
 
 // forward
 class particles;
-class dictionary;
 
 /**
  * Base class for particle insertion
  */
-class insertion
+class insertion : public fileDictionary
 {
-protected:
+private:
 
 	/// Is insertion active
 	Logical    active_ = "No";
@@ -44,14 +41,20 @@ protected:
 	/// Check for collision? It is not active now
 	Logical    checkForCollision_ = "No";
 
+	/// if increase velocity in case particles are failed
+	/// to be inserted due to collision
+	Logical    increaseVelocity_ = "No";
+
 	/// Ref to particles
 	particles& particles_;
 
+	bool       readFromFile_ = false;
+
 	/// Read from dictionary
 	bool       readInsertionDict(const dictionary& dict);
 
 	/// Write to dictionary
-	bool writeInsertionDict(dictionary& dict)const;
+	// bool writeInsertionDict(dictionary& dict)const;
 
 public:
 
@@ -59,24 +62,47 @@ public:
 	TypeInfo("insertion");
 
 	/// Construct from component
-	insertion(particles& prtcl);
+	explicit insertion(particles& prtcl);
 
 	/// Destructor
-	virtual ~insertion() = default;
+	~insertion() override = default;
 
 	/// is Insertion active
-	bool isActive()const {
+	inline bool isActive() const
+	{
 		return active_();
 	}
 
-	/// read from iIstream
+	inline bool checkForCollision() const
+	{
+		return checkForCollision_();
+	}
+
+	inline particles& Particles()
+	{
+		return particles_;
+	}
+
+	inline const particles& Particles() const
+	{
+		return particles_;
+	}
+
+	inline bool readFromFile() const
+	{
+		return readFromFile_;
+	}
+
+	/// read from stream
+	bool read(iIstream& is, const IOPattern& iop) override;
+
+	/*/// read from iIstream
 	virtual bool read(iIstream& is) = 0;
 
 	/// write to iOstream
-	virtual bool write(iOstream& os)const = 0;
+	virtual bool write(iOstream& os)const = 0;*/
-
-
 };
+
 }
 
 #endif

src/Particles/Insertion/insertionRegion/insertionRegion.cpp

@@ -18,42 +18,79 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #include "insertionRegion.hpp"
 #include "dictionary.hpp"
+#include "insertion.hpp"
+#include "particles.hpp"
+#include "twoPartEntry.hpp"
+#include "types.hpp"
 
-bool pFlow::insertionRegion::readInsertionRegion
+namespace pFlow
-(
-	const dictionary& dict
-)
 {
+template<typename T>
+bool
+setOneEntry(const twoPartEntry& tpEntry, anyList& varList)
+{
+	if (getTypeName<T>() != tpEntry.firstPart())
+		return false;
 
-	name_ = dict.name();
+	T val = tpEntry.secondPartVal<T>();
-	type_ = dict.getVal<word>("type");
+	varList.emplaceBack(tpEntry.keyword(), val);
+
+	return true;
+}
+
+bool
+readOneEtrty(const dataEntry& entry, anyList& varList)
+{
+	twoPartEntry stField(entry);
+
+	if (!(setOneEntry<real>(stField, varList) ||
+	      setOneEntry<realx3>(stField, varList) ||
+	      setOneEntry<realx4>(stField, varList) ||
+	      setOneEntry<int8>(stField, varList) ||
+	      setOneEntry<uint8>(stField, varList) ||
+	      setOneEntry<uint32>(stField, varList) ||
+	      setOneEntry<uint64>(stField, varList) ||
+	      setOneEntry<int32>(stField, varList) ||
+	      setOneEntry<int64>(stField, varList)))
+	{
+		fatalErrorInFunction << "un-supported data type " << stField.firstPart()
+		                     << endl;
+		return false;
+	}
+	return true;
+}
+}
+
+bool
+pFlow::insertionRegion::readInsertionRegion(const dictionary& dict)
+{
+	type_ = dict.getVal<word>("regionType");
+
+	rate_ = dict.getVal<real>("rate");
 
 	pRegion_ = peakableRegion::create(type_, dict.subDict(type_ + "Info"));
 
-	mixture_ = makeUnique<shapeMixture>(dict.subDict("mixture"));
+	mixture_ = makeUnique<shapeMixture>(
+	  dict.subDict("mixture"),
+	  insertion_.Particles().getShapes().shapeNameList()
+	);
 
-	addToNumInserted(mixture_().totalInserted());
+	numInserted_ = mixture_().totalInserted();
 
-	if( !dict.containsDictionay("setFields"))
+	if (dict.containsDictionay("setFields"))
 	{
-		output<<"\n insertion region "<< name_ << " does not contain setFields dictionary."
+		setFieldDict_ =
-		" An empty dictoiinary is created for it. \n";
+		  makeUnique<dictionary>("setFields", dict, dict.subDict("setFields"));
-		setFields_ = makeUnique<setFieldList>( dictionary("setFields") );
-	}
-	else
-	{
-		setFields_ = makeUnique<setFieldList>( dict.subDict("setFields") );
 	}
 
-	for(auto& sfEntry:setFields_())
+	if (setFieldDict_)
 	{
-		if(!sfEntry.checkForTypeAndValueAll())
+		if (!readSetFieldDict())
 		{
-			fatalErrorInFunction<<
+			fatalErrorInFunction << "Error in reading dictionary "
-			"   error in setFields dictionary "<< dict.globalName()<<endl;
+			                     << setFieldDict_().globalName() << endl;
 			return false;
 		}
 	}
@@ -61,79 +98,84 @@ bool pFlow::insertionRegion::readInsertionRegion
 	return true;
 }
 
-bool pFlow::insertionRegion::writeInsertionRegion
+bool
-(
+pFlow::insertionRegion::writeInsertionRegion(dictionary& dict) const
-	dictionary& dict
-) const
 {
-	
+	if (!dict.add("type", type_))
-	if(!dict.add("type", type_)) return false;
+		return false;
-	
 
 	if (pRegion_)
 	{
 		auto& prDict = dict.subDictOrCreate(type_ + "Info");
-		if(!pRegion_().write(prDict)) return false;
+		if (!pRegion_().write(prDict))
+			return false;
 	}
 
 	if (mixture_)
 	{
 		auto& mixDict = dict.subDictOrCreate("mixture");
-		if(!mixture_().write(mixDict)) return false;
+		if (!mixture_().write(mixDict))
+			return false;
 	}
 
-	if(setFields_)
+	/*if(setFields_)
 	{
 	    auto& sfDict = dict.subDictOrCreate("setFields");
 	    setFields_().write(sfDict);
+	}*/
+
+	return true;
+}
+
+bool
+pFlow::insertionRegion::readSetFieldDict()
+{
+	wordList Keys = setFieldDict_().dataEntryKeywords();
+
+	for (const auto& key : Keys)
+	{
+		if (!readOneEtrty(setFieldDict_().dataEntryRef(key), setFieldList_))
+		{
+			return false;
+		}
 	}
 
 	return true;
 }
 
-pFlow::insertionRegion::insertionRegion
+pFlow::insertionRegion::insertionRegion(
-(
+  const word&      name,
-	const dictionary& dict
+  const insertion& instn
 )
-:
+  : name_(name),
-	timeFlowControl(dict)
+    dict_(instn.subDict(name)),
+    insertion_(instn),
+    tControl_(dict_, "insertion")
 {
-
+	if (!readInsertionRegion(dict_))
-	if(!readInsertionRegion(dict))
 	{
 		fatalExit;
 	}
 }
 
-pFlow::insertionRegion::insertionRegion
+pFlow::uint32
-(
+pFlow::insertionRegion::numberToBeInserted(uint32 iter, real t, real dt)
-	const insertionRegion& src
-)
-:
-	timeFlowControl(src),
-	name_(src.name_),
-	type_(src.type_),
-	pRegion_( src.pRegion_? src.pRegion_->clone(): nullptr),
-	mixture_( src.mixture_? src.mixture_->clone(): nullptr),
-	setFields_( src.setFields_? src.setFields_->clone(): nullptr)
-{}
-
-pFlow::insertionRegion& pFlow::insertionRegion::operator=
-(
-	const insertionRegion& src
-)
 {
+	if (!tControl_.isInRange(iter, t, dt))
+		return 0u;
 
-	if(&src == this)return *this;
+	if (tControl_.isTimeStep())
-	timeFlowControl::operator=(src);
+	{
-
+		return static_cast<uint32>(
-	name_ = src.name_;
+		  (iter - tControl_.startIter() + tControl_.iInterval()) * dt * rate_ -
-	type_ = src.type_;
+		  numInserted_
-	pRegion_ = (src.pRegion_? src.pRegion_->clone(): nullptr);
+		);
-	mixture_ = (src.mixture_? src.mixture_->clone(): nullptr);
+	}
-	setFields_ = (src.setFields_? src.setFields_->clone(): nullptr);
+	else
-
+	{
-	return *this;
+		return static_cast<uint32>(
+		  (t - tControl_.startTime() + tControl_.rInterval()) * rate_ -
+		  numInserted_
+		);
+	}
 }
-
-

src/Particles/Insertion/insertionRegion/insertionRegion.hpp

@@ -21,15 +21,16 @@ Licence:
 #ifndef __insertionRegion_hpp__
 #define __insertionRegion_hpp__
 
-#include "timeFlowControl.hpp"
+#include "anyList.hpp"
+#include "baseTimeControl.hpp"
+#include "peakableRegion.hpp"
 #include "shapeMixture.hpp"
-#include "peakableRegions.hpp"
-#include "setFieldList.hpp"
 
 namespace pFlow
 {
 
 class dictionary;
+class insertion;
 
 /**
  * This class defines all the necessary enteties for defining an insertion
@@ -54,7 +55,8 @@ class dictionary;
 
     setFields
     {
-		velocity    realx3 (0.0 0.0 -0.6); // initial velocity of inserted particles 
+        velocity    realx3 (0.0 0.0 -0.6); // initial velocity of inserted
+particles
     }
 
     mixture
@@ -75,18 +77,32 @@ class dictionary;
  * | endTime | real | end of insertion (s) | No  |
  * | interval | real | time interval between successive insertions (s) | No |
  * | ###Info | dictionary | data for insertion region | No |
- * | setFields | dictionary | set field for inserted particles (s) | Yes [empty dictionray] |
+ * | setFields | dictionary | set field for inserted particles (s) | Yes [empty
+dictionray] |
  * | mixture | dictionary | mixture of particles to be inserted (s) | No |
  *
  */
 class insertionRegion
-:
-	public timeFlowControl
 {
-protected:
+private:
 
-	/// name of the region 
+	/// name of this region
-	word  		name_;
+	const word                name_;
+
+	/// insertion region dictionary
+	const dictionary&         dict_;
+
+	/// ref to pointStructure
+	const insertion&          insertion_;
+
+	/// @brief time control for insertion events
+	baseTimeControl           tControl_;
+
+	/// rate of insertion
+	real                      rate_;
+
+	/// number of inserted particles
+	uint32                    numInserted_ = 0;
 
 	/// type of insertion region
 	word                      type_;
@@ -97,9 +113,15 @@ protected:
 	/// mixture of shapes
 	uniquePtr<shapeMixture>   mixture_ = nullptr;
 
-	/// setFields for insertion region 
+	/// @brief  dictionary for set field
-	uniquePtr<setFieldList>     setFields_ = nullptr;
+	uniquePtr<dictionary>     setFieldDict_ = nullptr;
 
+	/// list of (filedName type value) for the fields
+	anyList                   setFieldList_;
+
+private:
+
+	// - private methods
 
 	/// read from dictionary
 	bool readInsertionRegion(const dictionary& dict);
@@ -107,6 +129,7 @@ protected:
 	/// write to dictionary
 	bool writeInsertionRegion(dictionary& dict) const;
 
+	bool readSetFieldDict();
 
 public:
 
@@ -116,44 +139,68 @@ public:
 	// - Constructors
 
 	/// Construct from a dictionary
-		insertionRegion(const dictionary& dict);
+	insertionRegion(const word& name, const insertion& instn);
-
-		/// Copy
-		insertionRegion(const insertionRegion& src);
-
-		/// Move 
-		insertionRegion(insertionRegion&&) = default;
-
-		/// Copy assignment
-		insertionRegion& operator=(const insertionRegion&);
-
-		/// Move assignment
-		insertionRegion& operator=(insertionRegion&&) = default;
 
 	/// Destructor
 	~insertionRegion() = default;
 
-
 	// - Methods
 
-		/// Const ref to setFields
-		const auto& setFields()const
-		{
-			return setFields_();
-		}
-
 	/// Const ref to name of the region
 	const auto& name() const
 	{
 		return name_;
 	}
 
+	/// return type of insertion region
+	const auto& type() const
+	{
+		return type_;
+	}
+
+	const auto& Insertion() const
+	{
+		return insertion_;
+	}
+
+	inline bool insertionTime(uint32 iter, real t, real dt) const
+	{
+		return tControl_.timeEvent(iter, t, dt);
+	}
+
+	uint32        numberToBeInserted(uint32 iter, real t, real dt);
+
+	inline uint32 addToNumInserted(uint32 newInserted)
+	{
+		return numInserted_ += newInserted;
+	}
+
+	inline uint32 totalInserted() const
+	{
+		return numInserted_;
+	}
+
+	auto& mixture()
+	{
+		return mixture_();
+	}
+
+	auto& pRegion()
+	{
+		return pRegion_();
+	}
+
+	const auto& setFieldList() const
+	{
+		return setFieldList_;
+	}
 	// - IO operation
 
 	/// read from dictionary
-		bool read(const dictionary& dict)
+	/*bool read(const dictionary& dict)
 	{
-			if(!timeFlowControl::read(dict))return false;
+	    if (!timeFlowControl::read(dict))
+	        return false;
 
 	    return readInsertionRegion(dict);
 	}
@@ -161,10 +208,11 @@ public:
 	/// write to dictionary
 	bool write(dictionary& dict) const
 	{
-			if(!timeFlowControl::write(dict)) return false;
+	    if (!timeFlowControl::write(dict))
+	        return false;
 
 	    return writeInsertionRegion(dict);
-		}
+	}*/
 };
 
 } // pFlow

src/Particles/Insertion/insertionRegion/timeFlowControl.cpp

@@ -1,83 +0,0 @@
-/*------------------------------- phasicFlow ---------------------------------
-      O        C enter of
-     O O       E ngineering and
-    O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
-------------------------------------------------------------------------------
-  Copyright (C): www.cemf.ir
-  email: hamid.r.norouzi AT gmail.com
-------------------------------------------------------------------------------  
-Licence:
-  This file is part of phasicFlow code. It is a free software for simulating 
-  granular and multiphase flows. You can redistribute it and/or modify it under
-  the terms of GNU General Public License v3 or any other later versions. 
- 
-  phasicFlow is distributed to help others in their research in the field of 
-  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
-  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-
------------------------------------------------------------------------------*/
-
-#include "timeFlowControl.hpp"
-#include "dictionary.hpp"
-
-size_t pFlow::timeFlowControl::numberToBeInserted(real currentTime)
-{
-	if(currentTime<startTime_)return 0;
-	if(currentTime>endTime_) return 0;
-
-	return static_cast<size_t>
-	( 
-		(currentTime - startTime_ + interval_)*rate_ - numInserted_ 
-	);
-}
-
-bool pFlow::timeFlowControl::readTimeFlowControl
-(
-	const dictionary& dict
-)
-{
-	rate_ = dict.getVal<real>("rate");
-	startTime_ 	= dict.getVal<real>("startTime");
-	endTime_	= dict.getVal<real>("endTime");
-	interval_	= dict.getVal<real>("interval");
-	numInserted_=0;
-	return true;
-}
-
-bool pFlow::timeFlowControl::writeTimeFlowControl
-(
-	dictionary& dict
-) const
-{
-	if(!dict.add("rate", rate_)) return false;
-	if(!dict.add("startTime", startTime_)) return false;
-	if(!dict.add("endTime", endTime_)) return false;
-	if(!dict.add("interval", interval_)) return false;
-
-	return true;
-}
-
-
-pFlow::timeFlowControl::timeFlowControl
-(
-	const dictionary& dict
-)
-{
-
-	if(!readTimeFlowControl(dict))
-	{
-		fatalExit;
-	}
-
-}
-
-bool pFlow::timeFlowControl::insertionTime( real currentTime, real dt)
-{
-	if(currentTime < startTime_) return false;
-
-	if(currentTime > endTime_) return false;
-	if( mod(abs(currentTime-startTime_),interval_)/dt < 1 ) return true;
-
-	return false;
-}

src/Particles/Insertion/insertionRegion/timeFlowControl.hpp

@@ -1,101 +0,0 @@
-/*------------------------------- phasicFlow ---------------------------------
-      O        C enter of
-     O O       E ngineering and
-    O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
-------------------------------------------------------------------------------
-  Copyright (C): www.cemf.ir
-  email: hamid.r.norouzi AT gmail.com
-------------------------------------------------------------------------------  
-Licence:
-  This file is part of phasicFlow code. It is a free software for simulating 
-  granular and multiphase flows. You can redistribute it and/or modify it under
-  the terms of GNU General Public License v3 or any other later versions. 
- 
-  phasicFlow is distributed to help others in their research in the field of 
-  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
-  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-
------------------------------------------------------------------------------*/
-
-#ifndef __timeFlowControl_hpp__ 
-#define __timeFlowControl_hpp__ 
-
-#include "types.hpp"
-#include "streams.hpp"
-
-namespace pFlow
-{
-
-class dictionary;
-
-/**
- *  Time control for particle insertion
- */
-class timeFlowControl
-{
-protected:
-
-	/// start time of insertion
-	real startTime_;
-
-	/// end time of insertion 
-	real endTime_;
-
-	/// time interval between each insertion
-	real interval_;
-
-	/// rate of insertion 
-	real rate_;
-
-	/// number of inserted particles
-	size_t numInserted_ = 0;	
-
-
-	/// Read dictionary
-	bool readTimeFlowControl(const dictionary& dict);
-
-	/// Write to dictionary 
-	bool writeTimeFlowControl(dictionary& dict) const;
-
-	/// Return number of particles to be inserted at time currentTime
-	size_t numberToBeInserted(real currentTime);
-	
-	/// Add to numInserted
-	inline
-	size_t addToNumInserted(size_t newInserted)
-	{
-		return numInserted_ += newInserted;
-	}
-
-public:
-
-	/// Construct from dictionary 
-	timeFlowControl(const dictionary& dict);
-
-	/// Is currentTime the insertion moment?
-	bool insertionTime( real currentTime, real dt);
-	
-	/// Total number inserted so far
-	size_t totalInserted()const
-	{
-		return numInserted_;
-	}
-
-	/// Read from dictionary
-	bool read(const dictionary& dict)
-	{
-		return readTimeFlowControl(dict);
-	}
-
-	/// Write to dictionary 
-	bool write(dictionary& dict)const
-	{
-		return writeTimeFlowControl(dict);
-	}
-
-};
-
-}
-
-#endif //__timeFlowControl_hpp__ 

src/Particles/Insertion/shapeMixture/shapeMixture.cpp

@@ -18,25 +18,44 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #include "shapeMixture.hpp"
 #include "dictionary.hpp"
 
-
+pFlow::shapeMixture::shapeMixture(
-pFlow::shapeMixture::shapeMixture
+  const dictionary& dict,
-(
+  const wordList&   validNames
-	const dictionary & dict
 )
 {
 	if (!read(dict))
 	{
 		fatalExit;
 	}
+
+	for (const auto& rN : names_)
+	{
+		bool found = false;
+		for (const auto& vN : validNames)
+		{
+			if (rN == vN)
+			{
+				found = true;
+				break;
+			}
 		}
 
-pFlow::word pFlow::shapeMixture::getNextShapeName()
+		if (!found)
 		{
+			fatalErrorInFunction
+			  << "Shape name " << rN << " provided in mixture dictionary "
+			  << dict.globalName() << " is invalid. \n Valid names are "
+			  << validNames << endl;
+		}
+	}
+}
 
+pFlow::word
+pFlow::shapeMixture::getNextShapeName()
+{
 	ForAll(i, names_)
 	{
 		if (current_[i] < number_[i])
@@ -51,24 +70,20 @@ pFlow::word pFlow::shapeMixture::getNextShapeName()
 	return getNextShapeName();
 }
 
-void pFlow::shapeMixture::getNextShapeNameN
+void
-(
+pFlow::shapeMixture::getNextShapeNameN(size_t n, wordVector& names)
-	size_t n,
-	wordVector& names
-)
 {
 	names.clear();
 
-
+	for (size_t i = 0u; i < n; ++i)
-	for(label i=0; i<n; ++i)
 	{
 		names.push_back(getNextShapeName());
 	}
 }
 
-bool pFlow::shapeMixture::read(const dictionary& dict)
+bool
+pFlow::shapeMixture::read(const dictionary& dict)
 {
-
 	bool containNumberIneserted = false;
 
 	auto shNames = dict.dataEntryKeywords();
@@ -90,9 +105,10 @@ bool pFlow::shapeMixture::read(const dictionary& dict)
 		uint32 num = dict.getVal<uint32>(nm);
 		if (num <= 0)
 		{
-    		fatalErrorInFunction<<
+			fatalErrorInFunction << "  number inserte in front of " << nm
-    		"  number inserte in front of "<< nm << " is invalid: "<< num<<endl<<
+			                     << " is invalid: " << num << endl
-    		"  in dictionary "<<dict.globalName()<<endl;
+			                     << "  in dictionary " << dict.globalName()
+			                     << endl;
 			return false;
 		}
 		number_.push_back(num);
@@ -104,14 +120,16 @@ bool pFlow::shapeMixture::read(const dictionary& dict)
 	}
 	else
 	{
-    	numberInserted_ = uint32Vector(size(), static_cast<uint32>(0));
+		numberInserted_ =
+		  uint32Vector(numberInserted_.name(), size(), static_cast<uint32>(0));
 	}
 
 	if (numberInserted_.size() != names_.size())
 	{
-    	fatalErrorInFunction<<
+		fatalErrorInFunction
-    	"  number of elements in numberInserted ("<<numberInserted_.size()<<
+		  << "  number of elements in numberInserted ("
-    	") is not equal to the number of shape names: "<< names_<<endl;
+		  << numberInserted_.size()
+		  << ") is not equal to the number of shape names: " << names_ << endl;
 		return false;
 	}
 
@@ -124,26 +142,25 @@ bool pFlow::shapeMixture::read(const dictionary& dict)
 	return true;
 }
 
-bool pFlow::shapeMixture::write
+bool
-(
+pFlow::shapeMixture::write(dictionary& dict) const
-	dictionary& dict
-) const
 {
-
 	ForAll(i, names_)
 	{
 		if (!dict.add(names_[i], number_[i]))
 		{
-			fatalErrorInFunction<<
+			fatalErrorInFunction << "  error in writing " << names_[i]
-			"  error in writing "<< names_[i] << " to dictionary "<<dict.globalName()<<endl;
+			                     << " to dictionary " << dict.globalName()
+			                     << endl;
 			return false;
 		}
 	}
 
 	if (!dict.add("numberInserted", numberInserted_))
 	{
-		fatalErrorInFunction<<
+		fatalErrorInFunction
-		"  error in writing numberInserted to dictionary "<< dict.globalName()<<endl;
+		  << "  error in writing numberInserted to dictionary "
+		  << dict.globalName() << endl;
 		return false;
 	}
 

src/Particles/Insertion/shapeMixture/shapeMixture.hpp

@@ -21,8 +21,8 @@ Licence:
 #ifndef __shapeMixture_hpp__
 #define __shapeMixture_hpp__
 
-
 #include "Vectors.hpp"
+#include "Lists.hpp"
 
 namespace pFlow
 {
@@ -46,19 +46,19 @@ class dictionary;
  */
 class shapeMixture
 {
-protected:
+private:
 
 	/// List of shape names 
-	wordVector 		names_;
+	wordVector 		names_{"shapeNames"};
 
 	/// Number composition 
-	uint32Vector   	number_;
+	uint32Vector   	number_{"number"};
 
 	/// Number of inserted particles of each shape 
-	uint32Vector   	numberInserted_;
+	uint32Vector   	numberInserted_{"numberInserted"};
 
 	/// Current number of inserted
-	uint32Vector      current_;
+	uint32Vector      current_{"currentInserted"};
 
 public:
 
@@ -68,7 +68,7 @@ public:
 	// - Constrcutors 
 		
 		/// Construct from dictionary
-		shapeMixture(const dictionary & dict);
+		shapeMixture(const dictionary & dict, const wordList& validNames);
 
 		/// Copy
 		shapeMixture(const shapeMixture&) = default;
@@ -88,13 +88,6 @@ public:
 			return makeUnique<shapeMixture>(*this);
 		}
 
-		/// Polymorphic copy
-		shapeMixture* clonePtr()const
-		{
-			return new shapeMixture(*this);
-		}
-
-
 		/// Destructor 
 		~shapeMixture() = default;
 

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -198,7 +198,7 @@ bool pFlow::sphereParticles::initializeParticles()
 
 }*/
 
-bool pFlow::sphereParticles::initInertia()
+bool pFlow::sphereParticles::initializeParticles()
 {
 	
 	using exeSpace = typename realPointField_D::execution_space;
@@ -255,6 +255,55 @@ bool pFlow::sphereParticles::initInertia()
 	return true;
 }
 
+bool
+pFlow::sphereParticles::getParticlesInfoFromShape(
+  const wordVector& shapeNames,
+  uint32Vector&     propIds,
+  realVector&       diams,
+  realVector&       m,
+  realVector&       Is,
+  uint32Vector&     shIndex
+)
+{
+	auto numNew = static_cast<uint32>(shapeNames.size());
+
+	propIds.clear();
+	propIds.reserve(numNew);
+
+	diams.clear();
+	diams.reserve(numNew);
+
+	m.clear();
+	m.reserve(numNew);
+
+	Is.clear();
+	Is.reserve(numNew);
+
+	shIndex.clear();
+	shIndex.reserve(numNew);
+	
+
+	for(const auto& name:shapeNames)
+	{
+		uint32 indx;
+		if(spheres_.shapeNameToIndex(name,indx))
+		{
+			shIndex.push_back(indx);
+			Is.push_back( spheres_.Inertia(indx));
+			m.push_back(spheres_.mass(indx));
+			diams.push_back(spheres_.boundingDiameter(indx));
+			propIds.push_back( spheres_.propertyId(indx));
+		}
+		else
+		{	
+			fatalErrorInFunction<<"Shape name "<< name << 
+			"does not exist. The list is "<<spheres_.shapeNameList()<<endl;
+			return false;
+		}
+	}
+	
+	return true;
+}
 
 pFlow::sphereParticles::sphereParticles(
 	systemControl &control,
@@ -364,7 +413,7 @@ pFlow::sphereParticles::sphereParticles(
 
 	WARNING<<"setFields for rVelIntegration_"<<END_WARNING;
 	
-	if(!initInertia())
+	if(!initializeParticles())
 	{
 		fatalErrorInFunction;
 		fatalExit;
@@ -504,6 +553,59 @@ bool pFlow::sphereParticles::iterate()
 	return true;
 }
 
+bool pFlow::sphereParticles::insertParticles
+(
+	const realx3Vector &position, 
+	const wordVector &shapesNames, 
+	const anyList &setVarList
+)
+{
+	anyList newVarList(setVarList);
+
+	realVector mass("mass");
+	realVector I("I");
+	realVector diameter("diameter");
+	uint32Vector propId("propId");
+	uint32Vector shapeIdx("shapeIdx");
+
+	if(!getParticlesInfoFromShape(
+		shapesNames,
+		propId,
+		diameter,
+		mass,
+		I,
+		shapeIdx))
+	{
+		return false;
+	}
+
+	newVarList.emplaceBack(
+		mass_.name()+"Vector",
+		std::move(mass));
+
+	newVarList.emplaceBack(
+		I_.name()+"Vector",
+		std::move(I));
+	
+	newVarList.emplaceBack(
+		diameter_.name()+"Vector",
+		std::move(diameter));
+
+	newVarList.emplaceBack(
+		propertyId_.name()+"Vector",
+		std::move(propId));
+
+	newVarList.emplaceBack(
+		shapeIndex().name()+"Vector",
+		std::move(shapeIdx));
+	
+	if(!dynPointStruct().insertPoints(position, newVarList))
+	{
+		return false;
+	}
+
+    return true;
+}
 
 pFlow::word pFlow::sphereParticles::shapeTypeName()const
 {

src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp

@@ -29,23 +29,21 @@ Licence:
 #ifndef __sphereParticles_hpp__
 #define __sphereParticles_hpp__
 
-#include "systemControl.hpp"
-#include "particles.hpp"
-#include "sphereShape.hpp"
-#include "property.hpp"
 #include "indexContainer.hpp"
-
+#include "particles.hpp"
-
+#include "property.hpp"
+#include "sphereShape.hpp"
+#include "systemControl.hpp"
 
 namespace pFlow
 {
 
-class sphereParticles
+class sphereParticles : public particles
-:
-	public particles 
 {
 public:
+
 	using ShapeType = sphereShape;
+
 private:
 
 	/// reference to shapes
@@ -81,8 +79,19 @@ private:
 	/// timer for integration computations (correction step)
 	Timer                  intCorrectTimer_;
 
-	bool initInertia();
 	
+
+private:
+	bool initializeParticles();
+
+	bool getParticlesInfoFromShape(
+	  const wordVector& shapeNames,
+	  uint32Vector&     propIds,
+	  realVector&       diams,
+	  realVector&       m,
+	  realVector&       Is,
+	  uint32Vector&     shIndex
+	);
 	/*bool initializeParticles();
 
 	bool insertSphereParticles(
@@ -96,19 +105,15 @@ private:
 public:
 
 	/// construct from systemControl and property
-	sphereParticles(
+	sphereParticles(systemControl& control, const property& prop);
-		systemControl &control, 
-		const property& prop);
-
 
 	~sphereParticles() override = default;
 
 	/**
 	 * Insert new particles in position with specified shapes
 	 *
-	 * This function is involked by inserted object to insert new set of particles
+	 * This function is involked by inserted object to insert new set of
-	 * into the simulation. 
+	 * particles into the simulation. \param position position of new particles
-	 * \param position position of new particles
 	 * \param shape shape of new particles
 	 * \param setField initial value of the selected fields for new particles
 	 */
@@ -147,8 +152,7 @@ public:
 		return rVelocity_;
 	}
 
-	bool hearChanges
+	bool hearChanges(
-	(
 	  real           t,
 	  real           dt,
 	  uint32         iter,
@@ -160,19 +164,16 @@ public:
 		return false;
 	}
 
-	
 	const uint32PointField_D& propertyId() const override
 	{
 		return propertyId_;
 	}
 
-	
 	const realPointField_D& diameter() const override
 	{
 		return diameter_;
 	}
 
-
 	const realPointField_D& mass() const override
 	{
 		return mass_;
@@ -184,6 +185,11 @@ public:
 	/// iterate particles
 	bool iterate() override;
 
+	bool insertParticles(
+	  const realx3Vector& position,
+	  const wordVector&   shapesNames,
+	  const anyList&      setVarList
+	) override;
 
 	realx3PointField_D& rAcceleration() override
 	{
@@ -204,7 +210,6 @@ public:
 
 	const shape& getShapes() const override;
 
-	
 	void boundingSphereMinMax(real& minDiam, real& maxDiam) const override;
 
 }; // sphereParticles

src/Particles/dynamicPointStructure/dynamicPointStructure.cpp

@@ -78,20 +78,6 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 
 }
 
-/*bool pFlow::dynamicPointStructure::beforeIteration()
-{
-	pointStructure::beforeIteration();
-	auto& acc = time().lookupObject<realx3PointField_D>("acceleration");
-    return predict(dt(), acc);
-}*/
-
-/*bool pFlow::dynamicPointStructure::iterate()
-{
-	pointStructure::iterate();
-	auto& acc = time().lookupObject<realx3PointField_D>("acceleration");
-    return correct(dt(), acc);
-    
-}*/
 
 bool pFlow::dynamicPointStructure::beforeIteration()
 {
@@ -137,74 +123,3 @@ bool pFlow::dynamicPointStructure::correct
 
 	return true;	
 }
-
-/*FUNCTION_H
-pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::dynamicPointStructure::insertPoints
-(
-	const realx3Vector& pos,
-	const List<eventObserver*>& exclusionList
-)
-{
-	auto newIndicesPtr = pointStructure::insertPoints(pos, exclusionList);
-
-	// no new point is inserted
-	if( !newIndicesPtr ) return newIndicesPtr;
-
-	if(!integrationPos_().needSetInitialVals()) return newIndicesPtr;
-
-	auto hVel = velocity_.hostView();
-	auto n = newIndicesPtr().size();
-	auto index = newIndicesPtr().indicesHost();
-
-	realx3Vector velVec(n, RESERVE());
-	for(auto i=0; i<n; i++)
-	{
-		velVec.push_back( hVel[ index(i) ]);
-	}
-
-	integrationPos_->setInitialVals(newIndicesPtr(), pos );
-	integrationVel_->setInitialVals(newIndicesPtr(), velVec );
-	
-	return newIndicesPtr;
-
-}*/
-
-
-/*bool pFlow::dynamicPointStructure::update(
-	const eventMessage& msg) 
-{
-	if( msg.isInsert())
-	{
-		
-		if(!integrationPos_->needSetInitialVals())return true;
-
-		const auto indexHD = pStruct().insertedPointIndex();
-
-	
-		auto n = indexHD.size();
-
-		if(n==0) return true;
-
-		auto index = indexHD.indicesHost();
-
-		realx3Vector pos(n,RESERVE());
-		realx3Vector vel(n,RESERVE());
-		const auto hVel = velocity().hostView();
-		const auto hPos = pStruct().pointPosition().hostView();
-
-		for(auto i=0; i<n; i++)
-		{
-			pos.push_back( hPos[index(i)]);
-			vel.push_back( hVel[index(i)]);
-		}
-
-
-		
-		integrationPos_->setInitialVals(indexHD, pos);
-
-		integrationVel_->setInitialVals(indexHD, vel);		
-
-	}
-
-	return true;
-}*/

src/Particles/dynamicPointStructure/dynamicPointStructure.hpp

@@ -87,11 +87,10 @@ public:
 	/// when the component should evolve along time.
 	bool iterate() override;
 	
-		
+	/// prediction step (if any), is called in beforeIteration	
-
-
 	bool predict(real dt, realx3PointField_D& acceleration);
 
+	/// correction step, is called in iterate 
 	bool correct(real dt, realx3PointField_D& acceleration);
 
 };

src/Particles/particles/particleIdHandler/particleIdHandler.cpp

@@ -20,25 +20,69 @@ Licence:
 
 #include "particleIdHandler.hpp"
 
-pFlow::particleIdHandler::particleIdHandler(uint32PointField_D &id)
+pFlow::particleIdHandler::particleIdHandler(pointStructure& pStruct)
 :
-	id_(id)
+	uint32PointField_D
+	(
+		objectFile
+		(
+			"id",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+		),
+		pStruct,
+		static_cast<uint32>(-1),
+		static_cast<uint32>(-1)
+	)
 {
+
+}
+
+bool
+pFlow::particleIdHandler::hearChanges(
+  real           t,
+  real           dt,
+  uint32         iter,
+  const message& msg,
+  const anyList& varList
+)
+{
+	if(msg.equivalentTo(message::ITEM_INSERT))
+	{
+		const word eventName = message::eventName(message::ITEM_INSERT);
+
+		const auto& indices = varList.getObject<uint32IndexContainer>(
+			eventName);
+
+		uint32 numNew = indices.size();
+		if(numNew == 0u)return true;
+
+		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
+	{
+		return uint32PointField_D::hearChanges(t,dt,iter, msg,varList);
+	}
 }
 
 pFlow::uniquePtr<pFlow::particleIdHandler>
-	pFlow::particleIdHandler::create(uint32PointField_D &id)
+pFlow::particleIdHandler::create(pointStructure& pStruct)
 {
     word idHType = angleBracketsNames(
         "particleIdHandler", 
         pFlowProcessors().localRunTypeName());
 
 
-	if( pointFieldvCtorSelector_.search(idHType) )
+	if( pointStructurevCtorSelector_.search(idHType) )
 	{
 		REPORT(1)<<"Creating particle id handler "<< 
 			Green_Text(idHType)<<END_REPORT;
-		return pointFieldvCtorSelector_[idHType] (id);
+		return pointStructurevCtorSelector_[idHType] (pStruct);
 	}
 	else
 	{
@@ -47,7 +91,7 @@ pFlow::uniquePtr<pFlow::particleIdHandler>
 			fatalError << "Ctor Selector "<< idHType << " dose not exist. \n"
 			<<"Avaiable ones are: \n\n"
 			,
-			pointFieldvCtorSelector_
+			pointStructurevCtorSelector_
 		);
 		fatalExit;
 	}

src/Particles/particles/particleIdHandler/particleIdHandler.hpp

@@ -28,34 +28,26 @@ namespace pFlow
 {
 
 class particleIdHandler
+:
+	public uint32PointField_D
 {
-private:
-	
-  uint32PointField_D&    id_;
-
-protected:
-
-	auto& id()
-	{
-		return id_;
-	}
   
 public:
 	
-	TypeInfo("particleIdHandler");
+	/// class info
+	ClassInfo("particleIdHandler");
 
-	explicit particleIdHandler(uint32PointField_D & id);
+	explicit particleIdHandler(pointStructure& pStruct);
 
 	create_vCtor
 	(
 		particleIdHandler,
-		pointField,
+		pointStructure,
-		(uint32PointField_D & id),
+		(pointStructure& pStruct),
-		(id)
+		(pStruct)
 	);
 
-	virtual 
+	~particleIdHandler()override=default;
-	~particleIdHandler()=default;
 
 	virtual
 	Pair<uint32, uint32> getIdRange(uint32 nNewParticles)=0;
@@ -63,9 +55,20 @@ public:
 	virtual 
 	bool initialIdCheck()=0;
 
+	// heat change for possible insertion of particles
+	// overrdie from internalField
+	bool hearChanges
+	(
+		real t,
+		real dt,
+		uint32 iter,
+		const message& msg, 
+		const anyList& varList
+	) override;
+
 	static
 	uniquePtr<particleIdHandler> create(
-		uint32PointField_D & id);
+		pointStructure& pStruct);
 	
 };
 

src/Particles/particles/particles.cpp

@@ -15,25 +15,16 @@ Licence:
   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 "particles.hpp"
 
-
+pFlow::particles::particles(systemControl& control)
-pFlow::particles::particles
+  : observer(defaultMessage_),
-(
-    systemControl &control
-)
-: 
-	observer(defaultMessage_),
     demComponent("particles", control),
     dynPointStruct_(control),
-	id_
+    /*id_(
-	(
+      objectFile(
-		objectFile
-		(
         "id",
         "",
         objectFile::READ_IF_PRESENT,
@@ -42,11 +33,9 @@ pFlow::particles::particles
       dynPointStruct_,
       static_cast<uint32>(-1),
       static_cast<uint32>(-1)
-	),
+    ),*/
-	shapeIndex_
+    shapeIndex_(
-	(
+      objectFile(
-		objectFile
-		(
         "shapeIndex",
         "",
         objectFile::READ_ALWAYS,
@@ -55,13 +44,13 @@ pFlow::particles::particles
       dynPointStruct_,
       0
     ),
-	accelertion_
+    accelertion_(
-	(
       objectFile(
         "accelertion",
         "",
         objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS),
+        objectFile::WRITE_ALWAYS
+      ),
       dynPointStruct_,
       zero3
     ),
@@ -70,51 +59,54 @@ pFlow::particles::particles
         "contactForce",
         "",
         objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS),
+        objectFile::WRITE_ALWAYS
+      ),
       dynPointStruct_,
-		zero3),
+      zero3
+    ),
     contactTorque_(
       objectFile(
         "contactTorque",
         "",
         objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS),
+        objectFile::WRITE_ALWAYS
+      ),
       dynPointStruct_,
-		zero3),
+      zero3
-	idHandler_(particleIdHandler::create(id_))
+    ),
+    idHandler_(particleIdHandler::create(dynPointStruct_))
 {
 	this->addToSubscriber(dynPointStruct_);
 
 	idHandler_().initialIdCheck();
-
 }
 
-bool pFlow::particles::beforeIteration() 
+bool
+pFlow::particles::beforeIteration()
 {
-	
 	zeroForce();
 	zeroTorque();
 
 	return dynPointStruct_.beforeIteration();
 }
 
-bool pFlow::particles::iterate()
+bool
+pFlow::particles::iterate()
 {
 	return dynPointStruct_.iterate();
 }
 
-bool pFlow::particles::afterIteration() 
+bool
+pFlow::particles::afterIteration()
 {
 	return dynPointStruct_.afterIteration();
 }
 
-void pFlow::particles::boundingSphereMinMax
+void
-(
+pFlow::particles::boundingSphereMinMax(real& minDiam, real& maxDiam) const
-	real &minDiam, 
-	real &maxDiam
-) const
 {
 	auto& shp = getShapes();
+
 	minDiam = shp.minBoundingSphere();
 	maxDiam = shp.maxBoundingSphere();
 }

src/Particles/particles/particles.hpp

@@ -1,49 +1,45 @@
-/*------------------------------- phasicFlow ---------------------------------
+/*------------------------------- phasicFlow
-      O        C enter of
+--------------------------------- O        C enter of O O       E
-     O O       E ngineering and
+ngineering and O   O      M ultiscale modeling of OOOOOOO     F luid
-    O   O      M ultiscale modeling of
+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
-  granular and multiphase flows. You can redistribute it and/or modify it under
+simulating granular and multiphase flows. You can redistribute it
-  the terms of GNU General Public License v3 or any other later versions. 
+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 
+  phasicFlow is distributed to help others in their research in the
-  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+field of granular and multiphase flows, but WITHOUT ANY WARRANTY;
-  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
 
-
 #ifndef __particles_hpp__
 #define __particles_hpp__
 
-#include "dynamicPointStructure.hpp"
 #include "demComponent.hpp"
-#include "shape.hpp"
+#include "dynamicPointStructure.hpp"
 #include "particleIdHandler.hpp"
+#include "shape.hpp"
 
 namespace pFlow
 {
 
 class particles
-:
+  : public observer
-	public observer,
+  , public demComponent
-	public demComponent
 {
 private:
 
 	/// dynamic point structure for particles center mass
 	dynamicPointStructure        dynPointStruct_;
 
-	/// id of particles on host
+	/// shape index of each particle 
-	uint32PointField_D		id_;
-
-	
 	uint32PointField_D           shapeIndex_;
 
 	/// acceleration on device
@@ -55,12 +51,13 @@ private:
 	/// contact torque field
 	realx3PointField_D           contactTorque_;
 
+	/// handling new ids for new particles 
 	uniquePtr<particleIdHandler> idHandler_ = nullptr;
 
-	static inline
+	/// messages for this objects
-	const message defaultMessage_{message::DEFAULT}; 
+	static inline const message  defaultMessage_{ message::DEFAULT };
-	
 
+protected:
 
 	void zeroForce()
 	{
@@ -72,8 +69,6 @@ private:
 		contactTorque_.fill(zero3);
 	}
 
-protected:
-
 	inline auto& dynPointStruct()
 	{
 		return dynPointStruct_;
@@ -84,14 +79,12 @@ protected:
 		return dynPointStruct_.pointPosition();
 	}
 
-	/*inline 
+	inline auto& idHandler()
-	auto& velocity()
 	{
-		return dynPointStruct_.velocity();
+		return idHandler_();
-	}*/
+	}
 	
-	inline 
+	inline auto& shapeIndex()
-	auto& shapeIndex()
 	{
 		return shapeIndex_;
 	}
@@ -103,42 +96,38 @@ public:
 
 	explicit particles(systemControl& control);
 
-
+	inline const auto& dynPointStruct() const
-	inline 
-	const auto& dynPointStruct()const 
 	{
 		return dynPointStruct_;
 	}
 
-	inline
+	inline const pointStructure& pStruct() const
-	const pointStructure& pStruct()const
 	{
 		return dynPointStruct_;
 	}
 
-	inline
+	inline const auto& simDomain() const
-	const auto& simDomain()const
 	{
 		return dynPointStruct_.simDomain();
 	}
 
-	inline 
+	inline const auto& thisDomain() const
-	const auto& thisDomain()const
 	{
 		return dynPointStruct_.thisDomain();
 	}
 
-	inline
+	inline const auto& extendedDomain() const
-	const auto& extendedDomain()const
 	{
 		return dynPointStruct_.extendedDomain();
 	}
 
-	inline auto size()const{
+	inline auto size() const
+	{
 		return dynPointStruct_.size();
 	}
 
-	inline auto capacity() const{
+	inline auto capacity() const
+	{
 		return dynPointStruct_.capacity();
 	}
 
@@ -152,50 +141,42 @@ public:
 		return dynPointStruct_.isAllActive();
 	}
 
-	inline 
+	inline const auto& pointPosition() const
-	const auto& pointPosition()const
 	{
 		return dynPointStruct_.pointPosition();
 	}
 
-	inline 
+	inline const auto& velocity() const
-	const auto& velocity()const
 	{
 		return dynPointStruct_.velocity();
 	}
 
-	inline 
+	inline const auto& accelertion() const
-	const auto& accelertion()const
 	{
 		return accelertion_;
 	}
 
-	inline 
+	inline auto& accelertion()
-	auto& accelertion()
 	{
 		return accelertion_;
 	}
 
-	inline 
+	inline auto& contactForce()
-	auto& contactForce() 
 	{
 		return contactForce_;
 	}
 
-	inline
+	inline const auto& contactForce() const
-	const auto& contactForce() const
 	{
 		return contactForce_;
 	}
 
-	inline
+	inline auto& contactTorque()
-	auto& contactTorque()
 	{
 		return contactTorque_;
 	}
 
-	inline
+	inline const auto& contactTorque() const
-	const auto& contactTorque() const 
 	{
 		return contactTorque_;
 	}
@@ -206,45 +187,32 @@ public:
 
 	bool         afterIteration() override;
 
-	/*virtual
+	virtual bool insertParticles(
-	bool insertParticles
-	(
 	  const realx3Vector& position,
-	 	const wordVector&  shapes,
+	  const wordVector&   shapesNames,
-	 	const setFieldList& setField
+	  const anyList&      setVarList
-	 ) = 0;*/
+	) = 0;
 
-	virtual 
+	virtual const uint32PointField_D& propertyId() const = 0;
-	const uint32PointField_D& propertyId()const = 0;
 
-	virtual
+	virtual const realPointField_D&   diameter() const = 0;
-	const realPointField_D& diameter()const = 0;
 
-	virtual 
+	virtual const realPointField_D&   mass() const = 0;
-	const realPointField_D& mass()const = 0;		
 
-	virtual
+	virtual realx3PointField_D&       rAcceleration() = 0;
-	realx3PointField_D& rAcceleration() = 0;
 
-	virtual
+	virtual const realx3PointField_D& rAcceleration() const = 0;
-	const realx3PointField_D& rAcceleration() const = 0;
 
-	virtual
+	virtual const realPointField_D&   boundingSphere() const = 0;
-	const realPointField_D& boundingSphere()const = 0;
 
-	virtual 
+	virtual word                      shapeTypeName() const = 0;
-	word shapeTypeName()const = 0;
-
-	virtual
-	const shape& getShapes()const = 0;
-
-	virtual
-	void boundingSphereMinMax(real & minDiam, real& maxDiam)const = 0;
 
+	virtual const shape&              getShapes() const = 0;
 
+	virtual void boundingSphereMinMax(real& minDiam, real& maxDiam) const = 0;
 
 }; // particles
 
-} // pFlow
+} // namespace pFlow
 
 #endif //__particles_hpp__

src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.cpp

@@ -1,18 +1,18 @@
-#include "regularParticleHandler.hpp"
+#include "regularParticleIdHandler.hpp"
 
 
 
-pFlow::regularParticleHandler::regularParticleHandler
+pFlow::regularParticleIdHandler::regularParticleIdHandler
 (
-    uint32PointField_D & id
+    pointStructure& pStruct
 )
 :
-    particleIdHandler(id)
+    particleIdHandler(pStruct)
 {
 }
 
 pFlow::Pair<pFlow::uint32, pFlow::uint32> 
-    pFlow::regularParticleHandler::getIdRange(uint32 nNewParticles)
+    pFlow::regularParticleIdHandler::getIdRange(uint32 nNewParticles)
 {
     uint32 startId;
     if(maxId_==-1)
@@ -28,18 +28,18 @@ pFlow::Pair<pFlow::uint32, pFlow::uint32>
     return {startId, endId};
 }
 
-bool pFlow::regularParticleHandler::initialIdCheck()
+bool pFlow::regularParticleIdHandler::initialIdCheck()
 {
     /// empty point structure / no particles in simulation
-    if( id().pStruct().empty() ) return true;
+    if( pStruct().empty() ) return true;
 
-    uint32 maxId = max(id());
+    uint32 maxId = max( *this );
     
     /// particles should get ids from 0 to size-1
     if(maxId == -1)
     {
-        fillSequence(id(),0u);
+        fillSequence(*this,0u);
-        maxId_ = id().size()-1;
+        maxId_ = size()-1;
     }
     else
     {

src/Particles/particles/regularParticleIdHandler/regularParticleIdHandler.hpp

@@ -17,15 +17,15 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
-#ifndef __regularParticleHandler_hpp__
+#ifndef __regularParticleIdHandler_hpp__
-#define __regularParticleHandler_hpp__
+#define __regularParticleIdHandler_hpp__
 
 #include "particleIdHandler.hpp"
 
 namespace pFlow
 {
 
-class regularParticleHandler
+class regularParticleIdHandler
 :
     public particleIdHandler
 {
@@ -35,17 +35,17 @@ private:
 
 public:
 
-    TypeInfo("particleIdHandler<regular>");
+    ClassInfo("particleIdHandler<regular>");
     
-    explicit regularParticleHandler(uint32PointField_D & id);
+    explicit regularParticleIdHandler(pointStructure& pStruct);
 
-    ~regularParticleHandler()override = default;
+    ~regularParticleIdHandler()override = default;
 
 	add_vCtor
 	(
 		particleIdHandler,
-		regularParticleHandler,
+		regularParticleIdHandler,
-		pointField
+		pointStructure
 	);
 
     Pair<uint32, uint32> getIdRange(uint32 nNewParticles)override;
@@ -57,4 +57,4 @@ public:
 }
 
 
-#endif //__regularParticleHandler_hpp__
+#endif //__regularParticleIdHandler_hpp__

src/Particles/particles/shape/shape.hpp

@@ -122,6 +122,17 @@ public:
 		return 0.0;
 	}
 
+	inline 
+	uint32 propertyId(uint32 index)const
+	{
+		if(indexValid(index))
+		{
+			return shapePropertyIds_[index];
+		}
+		fatalExit;
+		return 0.0;
+	}
+
 	virtual 
     real maxBoundingSphere()const = 0;
 

src/phasicFlow/CMakeLists.txt

@@ -5,6 +5,7 @@ types/basicTypes/Logical.cpp
 types/types.cpp
 
 globals/error.cpp
+globals/globalSettings.cpp
 
 processors/processors.cpp
 processors/localProcessors.cpp
@@ -62,11 +63,18 @@ containers/Field/Fields.cpp
 containers/symArrayHD/symArrays.cpp
 containers/List/anyList/anyList.cpp
 
+structuredData/box/box.cpp
+structuredData/sphere/sphere.cpp
+structuredData/cylinder/cylinder.cpp
+structuredData/peakableRegion/geometricRegion/geometricRegions.cpp
+structuredData/peakableRegion/PeakableRegion/PeakableRegions.cpp
+structuredData/peakableRegion/peakableRegion/peakableRegion.cpp
+
 structuredData/pointStructure/internalPoints/internalPointsKernels.cpp
 structuredData/pointStructure/internalPoints/internalPoints.cpp
 
 structuredData/zAxis/zAxis.cpp
-structuredData/box/box.cpp
+
 structuredData/line/line.cpp
 structuredData/infinitePlane/infinitePlane.cpp
 structuredData/plane/plane.cpp
@@ -82,10 +90,9 @@ structuredData/boundaries/boundaryReflective/boundaryReflective.cpp
 structuredData/boundaries/boundaryList.cpp
 structuredData/pointStructure/pointStructure/pointStructure.cpp
 structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp
-structuredData/pointStructure/selectors/selectBox/selectBox.cpp
 structuredData/pointStructure/selectors/selectRange/selectRange.cpp 
 structuredData/pointStructure/selectors/selectRandom/selectRandom.cpp
-
+structuredData/pointStructure/selectors/selectBox/selectBox.cpp
 
 triSurface/subSurface.cpp
 triSurface/triSurface.cpp

src/phasicFlow/Kokkos/Range.hpp

@@ -44,13 +44,11 @@ struct Range : public Kokkos::pair<T, T>
 	  //// - Constructors
 
 	  /// Default
-	  INLINE_FUNCTION_HD Range()
+	  INLINE_FUNCTION_HD Range() = default;
-	{
-	}
 
 	/// From end, set start to 0
 	INLINE_FUNCTION_HD
-	Range(const T& e)
+	explicit Range(const T& e)
 	  : Range(0, e)
 	{
 	}

src/phasicFlow/Kokkos/ViewAlgorithms.hpp

@@ -253,11 +253,14 @@ copy(
   uint32                                   sEnd
 )
 {
-	range32 sSpan(sStart, sEnd);
-	range32 dSpan(dStart, dStart + (sEnd - sStart));
 	
-	auto    srcSub = Kokkos::subview(src, sSpan);
+	auto srcSub = Kokkos::subview(
-	auto    dstSub = Kokkos::subview(dst, dSpan);
+		src, 
+		Kokkos::make_pair(sStart, sEnd));
+
+	auto dstSub = Kokkos::subview(
+		dst, 
+		Kokkos::make_pair(dStart, dStart + (sEnd - sStart)));
 
 	Kokkos::deep_copy(dstSub, srcSub);
 }

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

@@ -23,6 +23,7 @@ Licence:
 
 #include <any>
 
+#include "typeInfo.hpp"
 #include "types.hpp"
 #include "List.hpp"
 
@@ -44,8 +45,7 @@ public:
 
 	using const_iterator= typename anyListType::const_iterator;
 
-
+private:
-protected:
 
 	/// Contains a list of variables with any type 
 	anyListType 	anyList_;
@@ -109,7 +109,7 @@ public:
 				fatalExit;
 			}
 			names_.push_back(name);
-			types_.push_back(getTypeName<T>());
+			types_.push_back(getTypeName(other));
 			return anyList_.emplace_back(std::in_place_type<T>, other);
 		}
 
@@ -125,14 +125,23 @@ public:
 				fatalExit;
 			}
 			names_.push_back(name);
-			types_.push_back(getTypeName<T>());
+			types_.push_back(getTypeName(other));
-			return anyList_.emplace_back(std::in_place_type<T>, std::move(other));
+			return anyList_.emplace_back(std::in_place_type<T>, std::forward<T>(other));
 		}
 
 		/// Get the reference to variable by index
 		template<typename T>
 		T& getObject(size_t i)
 		{
+			if(getTypeName<T>() != types_[i])
+			{
+				fatalErrorInFunction<<
+				"requested object type is "<<getTypeName<T>()<<
+				" while strored type is "<<types_[i]<<
+				" for object "<< names_[i]<<endl;
+				fatalExit;
+			}
+
 			std::any *a = &(*anyList_.pos(i));
 			if(!a->has_value())
 			{
@@ -156,15 +165,7 @@ public:
 				"list of variables is "<<names_<<endl;
 				fatalExit;
 			}
-			std::any *a = &(*anyList_.pos(i));
+			return getObject<T>(static_cast<size_t>(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 
@@ -179,21 +180,22 @@ public:
 				"list of variables is "<<names_<<endl;
 				fatalExit;
 			}
-			const std::any *a = &(*anyList_.pos(i));
+			return getObject<T>(static_cast<size_t>(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
 		{
+			if(getTypeName<T>() != types_[i])
+			{
+				fatalErrorInFunction<<
+				"requested object type is "<<getTypeName<T>()<<
+				" while strored type is "<<types_[i]<<
+				" for object "<< names_[i]<<endl;
+				fatalExit;
+			}
+
 			const std::any *a = &(*anyList_.pos(i));
 			if(!a->has_value())
 			{
@@ -211,6 +213,15 @@ public:
 		{
 			if(i>= size())return nullptr;
 
+			if(getTypeName<T>() != types_[i])
+			{
+				warningInFunction<<
+				"requested object type is "<<getTypeName<T>()<<
+				" while strored type is "<<types_[i]<<
+				" for object "<< names_[i]<<endl;
+				return nullptr;
+			}
+
 			std::any *a = &(*anyList_.pos(i));
 			if( a->has_value())
 			{
@@ -227,6 +238,15 @@ public:
 		const T* getObjectPtr(size_t i)const
 		{
 			if(i>=size())return nullptr;
+			if(getTypeName<T>() != types_[i])
+			{
+				warningInFunction<<
+				"requested object type is "<<getTypeName<T>()<<
+				" while strored type is "<<types_[i]<<
+				" for object "<< names_[i]<<endl;
+				return nullptr;
+			}
+
 			const std::any *a = &(*anyList_.pos(i));
 			if( a->has_value())
 			{

src/phasicFlow/containers/VectorHD/VectorSingle.cpp

@@ -1,4 +1,3 @@
-#include "VectorSingle.hpp"
 /*------------------------------- phasicFlow ---------------------------------
       O        C enter of
      O O       E ngineering and
@@ -75,9 +74,7 @@ void pFlow::VectorSingle<T,MemorySpace>::changeCapacity
         {
             newView(i) = view_(i);
         }
-
         view_ = newView;
-
     }
     else
     {
@@ -661,7 +658,7 @@ auto pFlow::VectorSingle<T,MemorySpace>::getSpan()const
 
 template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H
-bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement(uint32IndexContainer indices, const T& val)
+bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement(const uint32IndexContainer& indices, const T& val)
 {
 	if(indices.empty()) return true;
 	
@@ -713,7 +710,7 @@ template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H
 bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
 (
-	const uint32IndexContainer indices, 
+	const uint32IndexContainer& indices, 
 	const std::vector<T>& vals
 )
 {
@@ -766,7 +763,7 @@ template<typename T, typename MemorySpace>
 INLINE_FUNCTION_H
 bool pFlow::VectorSingle<T,MemorySpace>::insertSetElement
 (
-	uint32IndexContainer indices, 
+	const uint32IndexContainer& indices, 
 	const ViewType1D<T, memory_space> vals
 )
 {

src/phasicFlow/containers/VectorHD/VectorSingle.hpp

@@ -90,10 +90,6 @@ private:
 
 	// - protected members and methods 		
 
-		/// growth factor for vector 
-		static const inline 
-		real 	growthFactor_ = vectorGrowthFactor__;
-
 		/// Is the memory of this vector accessible from Host
 		static constexpr 
 		bool  isHostAccessible_ = isHostAccessible<execution_space>();
@@ -108,7 +104,7 @@ private:
 	  	/// Evaluate capacity based on the input size 
 		static INLINE_FUNCTION_H uint32 evalCapacity(uint32 n)
 		{
-			return static_cast<uint32>(n*growthFactor_+1);
+			return static_cast<uint32>(n*pFlow::gSettings::vectorGrowthFactor__+1);
 		}
 
 	/// @brief Change size to n and preserve the conetent if realloc 
@@ -289,11 +285,6 @@ public:
 		INLINE_FUNCTION_H
 		void assign(const VectorType& src, bool srcCapacity = true);
 
-		
-		/*template<typename MSpace, 
-				std::enable_if_t<
-					!std::is_same_v<typename VectorSingle<T, MSpace>::memory_space, memory_space>, 
-					bool> = true>*/
 		template<typename MSpace>
 		INLINE_FUNCTION_H
 		void assignFromDevice(const VectorSingle<T, MSpace>& src, bool srcCapacity = true)
@@ -329,14 +320,14 @@ public:
 		auto getSpan()const;
 		
 		INLINE_FUNCTION_H
-		bool insertSetElement(uint32IndexContainer indices, const T& val);
+		bool insertSetElement(const uint32IndexContainer& indices, const T& val);
 		
 		INLINE_FUNCTION_H
-		bool insertSetElement(uint32IndexContainer indices, const std::vector<T>& vals);
+		bool insertSetElement(const uint32IndexContainer& indices, const std::vector<T>& vals);
 		
 		INLINE_FUNCTION_H
 		bool insertSetElement(
-			uint32IndexContainer indices, 
+			const uint32IndexContainer& indices, 
 			const ViewType1D<T, memory_space> vals);
 
 		INLINE_FUNCTION_H

src/phasicFlow/containers/indexContainer/indexContainer.hpp

@@ -24,7 +24,7 @@ Licence:
 #include <vector>
 
 #include "phasicFlowKokkos.hpp"
-
+#include "typeInfo.hpp"
 
 namespace pFlow
 {
@@ -94,6 +94,8 @@ protected:
 
 public:
 
+	TypeInfoTemplateNV11("indexContainer", IndexType);
+
 	//// - Constructors 
 
 		/// Default

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

@@ -17,11 +17,51 @@ Licence:
   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 -----------------------------------------------------------------------------*/
 
+
+template<class T, class MemorySpace>
+bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
+{
+	const word eventName = message::eventName(message::ITEM_INSERT);
+
+	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
+		T val = varList.getObject<T>(name());
+		success = field_.insertSetElement(indices, val);
+		
+	}
+	else if(varList.contains(name()+"Vector"))
+	{
+		// a vector of values is going to be assigned
+		const auto& valVec = varList.getObject<Vector<T>>(name()+"Vector");
+		success = field_.insertSetElement(indices,valVec);
+	}
+	else
+	{
+		success = field_.insertSetElement(indices, defaultValue_);
+	}
+
+	if(!success)
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	return true;
+
+}
+
 template<class T, class MemorySpace>
 pFlow::internalField<T, MemorySpace>::internalField
 (
 	const word& name, 
-	const internalPoints& internal
+	const internalPoints& internal,
+	const T& defVal
 )
 :
     observer
@@ -37,6 +77,7 @@ pFlow::internalField<T, MemorySpace>::internalField
 		internal.size(),
 		RESERVE()
 	),
+	defaultValue_(defVal),
 	internalPoints_(internal)
 {}
 
@@ -45,6 +86,7 @@ pFlow::internalField<T, MemorySpace>::internalField
 (
 	const word &name, 
 	const internalPoints &internal, 
+	const T& defVal,
 	const T &val
 )
 :
@@ -61,6 +103,7 @@ pFlow::internalField<T, MemorySpace>::internalField
 		internal.size(),
 		RESERVE()
 	),
+	defaultValue_(defVal),
 	internalPoints_(internal)
 {
 	fillInternal(val);
@@ -131,8 +174,7 @@ bool pFlow::internalField<T, MemorySpace>:: hearChanges
 	}
 	if(msg.equivalentTo(message::ITEM_INSERT))
 	{
-		notImplementedFunction;
+		return insert(varList);
-		return false;
 	}
 	return true;
 }

src/phasicFlow/containers/pointField/internalField/internalField.hpp

@@ -52,8 +52,13 @@ public:
 
 protected:
 
+	/// Field
     FieldType           	field_;
 
+	/// @brief value when a new item is added to field
+	T 						defaultValue_;
+
+	/// const ref to internal points 
     const internalPoints&   internalPoints_;
 
 	static inline
@@ -66,15 +71,19 @@ protected:
 		message::ITEM_DELETE
 	);
 
+	bool insert(const anyList& varList);
+
 public:
 
 	internalField(
 		const word& name, 
-		const internalPoints& internal);
+		const internalPoints& internal,
+		const T& defVal);
 	
 	internalField(
 		const word& name, 
 		const internalPoints& internal,
+		const T& defVal,
 		const T& val);
 	
 	inline 
@@ -142,6 +151,12 @@ public:
 		return field_.fieldKey();
 	}
 
+	inline
+	const T& defaultValue()const
+	{
+		return defaultValue_;
+	}
+
 	inline
 	auto activeRange()const
 	{

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

@@ -34,11 +34,11 @@ pFlow::pointField<T, MemorySpace>::pointField
 	InternalFieldType
 	(
 		objf.name(),
-		pStruct
+		pStruct,
+        defVal
 	),
 	boundaryFieldList_(pStruct.boundaries(), *this),
-	pStruct_(pStruct),
+	pStruct_(pStruct)
-	defaultValue_(defVal)
 {
     
     if(IOobject::implyRead())
@@ -77,11 +77,11 @@ pFlow::pointField<T, MemorySpace>::pointField
 	InternalFieldType
 	(
 		objf.name(),
-		pStruct
+		pStruct,
+        defVal
 	),
 	boundaryFieldList_(pStruct.boundaries(), *this),
-	pStruct_(pStruct),
+	pStruct_(pStruct)
-	defaultValue_(defVal)
 {
     if(IOobject::implyRead())
     {
@@ -125,11 +125,11 @@ pFlow::pointField<T, MemorySpace>::pointField
 	(
 		objf.name(),
 		pStruct,
+        defVal,
 		val
 	),
 	boundaryFieldList_(pStruct.boundaries(), *this),
-	pStruct_(pStruct),
+	pStruct_(pStruct)
-	defaultValue_(defVal)
 {
 
 }

src/phasicFlow/containers/pointField/pointField/pointField.hpp

@@ -61,9 +61,6 @@ private:
 		/// @brief  refrence to point structure 
 		const pointStructure&       pStruct_;
 
-		/// @brief value when a new item is added to field
-		T 	defaultValue_;
-
 
 public:
 	
@@ -96,6 +93,11 @@ public:
 
 	//// - Methods
 
+		const word& name() const final
+		{
+			return  IOobject::name();
+		}
+	
 		const auto& internal()const
 		{
 			return static_cast<const InternalFieldType&>(*this);

src/phasicFlow/eventManagement/message.hpp

@@ -108,6 +108,12 @@ public:
 		return *this;
 	}
 
+	inline const word& addAndName(EVENT evnt)
+	{
+		add(evnt);
+		return eventName(static_cast<size_t>(evnt));
+	}
+
 	inline message& remove(EVENT evnt)
 	{
 		events_.set(static_cast<size_t>(evnt), false);

src/phasicFlow/eventManagement/subscriber.hpp

@@ -76,6 +76,10 @@ public:
 		const message msg, 
 		const anyList& varList);
 	
+	const word& subscriberName()const
+	{
+		return subName_;
+	}
 
 };
 

src/phasicFlow/globals/globalSettings.cpp

@@ -0,0 +1,6 @@
+
+
+#include "globalSettings.hpp"
+
+
+const double pFlow::gSettings::vectorGrowthFactor__ = 1.1;

src/phasicFlow/globals/globalSettings.hpp

@@ -23,9 +23,11 @@ Licence:
 
 
 
-namespace pFlow
+namespace pFlow::gSettings
 {
-const inline double vectorGrowthFactor__ = 1.2;
+
+extern const double vectorGrowthFactor__;
+
 }
 
 

src/phasicFlow/globals/vocabs.hpp

@@ -25,31 +25,31 @@ namespace pFlow
 {
 
 // folders / repositories
-const inline char* settingsFolder__        = "settings";
+const inline char* const settingsFolder__        = "settings";
-const inline char* settingsRepository__    = "settings";
+const inline char* const settingsRepository__    = "settings";
-const inline char* caseSetupFolder__       = "caseSetup";
+const inline char* const caseSetupFolder__       = "caseSetup";
-const inline char* caseSetupRepository__   = "caseSetup";
+const inline char* const caseSetupRepository__   = "caseSetup";
-const inline char* geometryFolder__        = "geometry";
+const inline char* const geometryFolder__        = "geometry";
-const inline char* geometryRepository_     = "geometry";
+const inline char* const geometryRepository_     = "geometry";
-const inline char* integrationRepository__ = "integration";
+const inline char* const integrationRepository__ = "integration";
-const inline char* integrationFolder__     = "integration";
+const inline char* const integrationFolder__     = "integration";
 
 // file names
-const inline char* settingsFile__       = "settingsDict";
+const inline char* const settingsFile__       = "settingsDict";
-const inline char* domainFile__         = "domainDict";
+const inline char* const domainFile__         = "domainDict";
-const inline char* insertionFile__      = "particleInsertion";
+const inline char* const insertionFile__      = "particleInsertion";
-const inline char* shapeFile__          = "shapes";
+const inline char* const shapeFile__          = "shapes";
-const inline char* pointStructureFile__ = "pStructure";
+const inline char* const pointStructureFile__ = "pStructure";
-const inline char* triSurfaceFile__     = "triSurface";
+const inline char* const triSurfaceFile__     = "triSurface";
-const inline char* createParticles__    = "createParticles";
+const inline char* const createParticles__    = "createParticles";
-const inline char* motionModelFile__    = "motionModel";
+const inline char* const motionModelFile__    = "motionModel";
-const inline char* contactSearchFile__  = "contactSearch";
+const inline char* const contactSearchFile__  = "contactSearch";
-const inline char* propertyFile__       = "interaction";
+const inline char* const propertyFile__       = "interaction";
-const inline char* interactionFile__    = "interaction";
+const inline char* const interactionFile__    = "interaction";
-const inline char* postprocessFile__    = "postprocessDict";
+const inline char* const postprocessFile__    = "postprocessDict";
 
-const inline char* uniform__    = "uniform";
+const inline char* const uniform__    = "uniform";
-const inline char* nonUniform__ = "nonUniform";
+const inline char* const nonUniform__ = "nonUniform";
 
 }
 

src/phasicFlow/repository/Time/Time.hpp

@@ -28,6 +28,7 @@ Licence:
 #include "repository.hpp"
 #include "timeControl.hpp"
 
+
 namespace pFlow
 {
 

src/phasicFlow/repository/Time/timeControl.hpp

@@ -25,7 +25,7 @@ Licence:
 #include "types.hpp"
 #include "ranges.hpp"
 #include "streams.hpp"
-
+#include "timeInfo.hpp"
 
 namespace pFlow
 {
@@ -36,7 +36,7 @@ class dictionary;
 
 class timeControl
 {
-protected:
+private:
 
 
 	//// - Data members
@@ -90,7 +90,7 @@ protected:
 	
 public:
 
-	timeControl(const dictionary& dict);
+	explicit timeControl(const dictionary& dict);
 
 	timeControl(
 		dictionary& dict,
@@ -99,9 +99,7 @@ public:
 		real saveInterval,
 		word startTimeName);
 	
-	virtual ~timeControl()
+	virtual ~timeControl() = default;
-	{}
-
 	
 	real dt()const
 	{
@@ -186,7 +184,11 @@ public:
 		return timePrecision_;
 	}
 
-
+	inline
+	timeInfo TimeInfo()const
+	{
+		return {static_cast<uint32>(currentIter_), currentTime_, dt_};
+	}
 
 
 };

src/phasicFlow/repository/Time/timeInfo.hpp

@@ -15,31 +15,49 @@ Licence:
   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 __timeInfo__hpp_
+#define __timeInfo__hpp_
 
-#ifndef __peakableRegions_hpp__
+#include "types.hpp"
-#define __peakableRegions_hpp__
-
-#include "PeakableRegion.hpp"
-#include "boxRegion.hpp"
-#include "sphereRegion.hpp"
 
 namespace pFlow
 {
 
+class timeInfo
+{
+private:
 
-typedef PeakableRegion<boxRegion> boxPeakableRegion;
+	uint32 currentIter_;
 
-typedef PeakableRegion<sphereRegion> spherePeakableRegion;
+	real   currentTime_;
 
-} // pFlow
+	real   dt_;
 
-#ifndef BUILD_SHARED_LIBS
+public:
-	#include "peakableRegionInstantiate.cpp"
-#endif
 
+	timeInfo(uint32 cIter, real cTime, real dt)
+	  : currentIter_(cIter),
+	    currentTime_(cTime),
+	    dt_(dt)
+	{
+	}
 
+	inline const real& t() const
+	{
+		return currentTime_;
+	}
+	inline const real& dt() const
+	{
+		return dt_;
+	}
+	inline const uint32& iter() const
+	{
+		return currentIter_;
+	}
+};
 
-#endif
+} // namespace pFlow
+
+#endif //__timeInfo__hpp_

src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp

@@ -69,40 +69,7 @@ pFlow::uniquePtr<pFlow::pointField_H<Type>>
 	return Ptr;
 }
 
-/*template <typename Type>
-void* pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
-(
-	repository& owner,
-	pointStructure& pStruct,
-	bool verbose
-)
-{
 
-	if( !checkForType<Type>() ) return nullptr;
-	
-	Type defValue = entry_.secondPartVal<Type>();
-	
-	if(verbose)
-		REPORT(2)<<"Creating pointField " << Green_Text(fieldName())<< " with default value " << cyanText(defValue)<<
-	             " in repository "<< owner.name() <<END_REPORT;
-
-	// by default we perform operations on host
-	auto& field = 
-	owner.emplaceObject<pointField<Vector,Type, vecAllocator<Type>>>
-	(
-		objectFile
-		(
-			fieldName(),
-			"",
-			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
-		),
-		pStruct,
-		defValue
-	);
-
-	return &field;
-}*/
 
 template <typename Type>
 bool pFlow::setFieldEntry::setPointFieldSelected
@@ -148,45 +115,3 @@ bool pFlow::setFieldEntry::setPointFieldSelected
 	return false;
 }
 
-/*template <typename Type>
-void* pFlow::setFieldEntry::setPointFieldSelectedStd
-(
-	repository& owner,
-	int32IndexContainer& selected,
-	bool verbose
-)
-{
-	
-	if( !checkForType<Type>() ) return nullptr;
-
-	
-	auto fName = fieldName();
-
-	if( !owner.lookupObjectName(fName) )
-	{
-		fatalErrorInFunction<<
-		"  Cannot find "<< fName << " in repository " << owner.name() << ". \n";
-		return nullptr;
-	}
-
-	
-	Type value = entry_.secondPartVal<Type>();
-
-	if(verbose)
-		REPORT(2)<< "Setting selected points of " << Green_Text(fName)
-		 		 << " to value " << cyanText(value) <<END_REPORT;
-	
-	
-	auto fieldTypeName = owner.lookupObjectTypeName(fName);
-
-	if( pointField<Vector, Type, vecAllocator<Type>>::TYPENAME() == fieldTypeName )
-	{
-		auto& field = owner.lookupObject<pointField<Vector,Type, vecAllocator<Type>>>(fName);
-		if(field.insertSetElement(selected, value))
-			return &field;
-		else
-			return nullptr;
-	}
-
-	return nullptr;
-}*/

src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.cpp

@@ -19,8 +19,8 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 
-template<typename RegionType>
+template<typename GeomType>
-pFlow::PeakableRegion<RegionType>::PeakableRegion
+pFlow::PeakableRegion<GeomType>::PeakableRegion
 (
 	const word& type, 
 	const dictionary& dict
@@ -32,8 +32,8 @@ pFlow::PeakableRegion<RegionType>::PeakableRegion
 	
 }
 
-template<typename RegionType>
+template<typename GeomType>
-bool  pFlow::PeakableRegion<RegionType>::isInside
+bool  pFlow::PeakableRegion<GeomType>::isInside
 (
 	const realx3& point
 )const
@@ -42,20 +42,20 @@ bool  pFlow::PeakableRegion<RegionType>::isInside
 }
 
 
-template<typename RegionType>
+template<typename GeomType>
-pFlow::realx3 pFlow::PeakableRegion<RegionType>::peek()const
+pFlow::realx3 pFlow::PeakableRegion<GeomType>::peek()const
 {
 	return region_.peek();
 }
 	
-template<typename RegionType>
+template<typename GeomType>
-bool pFlow::PeakableRegion<RegionType>::read(const dictionary& dict)
+bool pFlow::PeakableRegion<GeomType>::read(const dictionary& dict)
 {
 	return region_.read(dict);
 }
 
-template<typename RegionType>
+template<typename GeomType>
-bool pFlow::PeakableRegion<RegionType>::write(dictionary& dict)const
+bool pFlow::PeakableRegion<GeomType>::write(dictionary& dict)const
 {
 	return region_.write(dict);
 }

src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegion.hpp

@@ -18,54 +18,49 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #ifndef __PeakableRegion_hpp__
 #define __PeakableRegion_hpp__
 
-
-#include "peakableRegion.hpp"
 #include "dictionary.hpp"
+#include "geometricRegion.hpp"
+#include "peakableRegion.hpp"
 
 namespace pFlow
 {
 
-
+template<typename GeomType>
-
+class PeakableRegion : public peakableRegion
-
-template<typename RegionType>
-class PeakableRegion
-:
-	public peakableRegion
 {
-protected:
+public:
+
+	using RegionType = geometricRegion<GeomType>;
+
+	using GeometryType = typename RegionType::GeometryType;
+
+private:
 
 	RegionType region_;
 
 public:
 
 	// type info
-	TypeInfoTemplate("peakableRegion", RegionType);
+	TypeInfoTemplate11("peakable", GeometryType);
 
 	//// - Constructors
 	PeakableRegion(const word& type, const dictionary& dict);
 
-		add_vCtor(
+	add_vCtor(peakableRegion, PeakableRegion, word);
-			peakableRegion,
-			PeakableRegion,
-			word
-			);
 
 	virtual uniquePtr<peakableRegion> clone() const override
 	{
-			return makeUnique<PeakableRegion<RegionType>>(*this);
+		return makeUnique<PeakableRegion<GeomType>>(*this);
 	}
 
 	virtual peakableRegion* clonePtr() const override
 	{
-			return new PeakableRegion<RegionType>(*this);
+		return new PeakableRegion<GeomType>(*this);
 	}
 
-
 	virtual ~PeakableRegion() = default;
 
 	//// - Methods
@@ -74,13 +69,11 @@ public:
 
 	virtual realx3 peek() const override;
 
-
 	//// - IO operatoins
 
 	virtual bool   read(const dictionary& dict) override;
 
 	virtual bool   write(dictionary& dict) const override;
-
 };
 
 } // pFlow

src/phasicFlow/structuredData/peakableRegion/PeakableRegion/PeakableRegions.cpp

@@ -19,9 +19,11 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 
-#include "peakableRegions.hpp"
+#include "PeakableRegion.hpp"
+#include "geometricRegions.hpp"
 
+template class pFlow::PeakableRegion<pFlow::box>;
 
-#ifdef BUILD_SHARED_LIBS
+template class pFlow::PeakableRegion<pFlow::sphere>;
-	#include "peakableRegionInstantiate.cpp"
+
-#endif
+template class pFlow::PeakableRegion<pFlow::cylinder>;

src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.cpp

@@ -0,0 +1,67 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+template<typename GeomType>
+pFlow::geometricRegion<GeomType>::geometricRegion(const dictionary& dict)
+  : geom_(dict)
+  , random_()
+{
+	minPoint_ = geom_.minPoint();
+	maxPoint_ = geom_.maxPoint();
+}
+
+template<typename GeomType>
+inline bool
+pFlow::geometricRegion<GeomType>::isInside(const realx3& p) const
+{
+	return geom_.isInside(p);
+}
+
+template<typename GeomType>
+inline pFlow::realx3
+pFlow::geometricRegion<GeomType>::peek() const
+{
+	for (uint32 i = 0u; i < maxTries_; i++)
+	{
+		realx3 p = random_.randomNumber(minPoint_, maxPoint_);
+
+		if (isInside(p))
+			return p;
+	}
+
+	fatalErrorInFunction << "cannot peek a random point from geometricRegion "
+	                    << typeName() << endl;
+	fatalExit;
+	return 0;
+}
+
+template<typename GeomType>
+bool
+pFlow::geometricRegion<GeomType>::read(const dictionary& dict)
+{
+	return geom_.read(dict);
+}
+
+template<typename GeomType>
+bool
+pFlow::geometricRegion<GeomType>::write(dictionary& dict) const
+{
+	return geom_.write(dict);
+}

src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegion.hpp

@@ -0,0 +1,75 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions.
+
+  phasicFlow is distributed to help others in their research in the field of
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __geometricRegion_hpp__
+#define __geometricRegion_hpp__
+
+
+#include "dictionary.hpp"
+#include "uniformRandomReal.hpp"
+#include "typeInfo.hpp"
+
+namespace pFlow
+{
+
+template<typename GeomType>
+class geometricRegion
+{
+public:
+
+	using GeometryType = GeomType;
+
+private:
+
+	GeomType                  geom_;
+
+	realx3                    minPoint_;
+
+	realx3                    maxPoint_;
+
+	mutable uniformRandomReal random_;
+
+	static constexpr uint32   maxTries_ = 100;
+
+public:
+
+	// - type info
+	TypeInfoTemplateNV11("geometricRegion", GeomType);
+
+	explicit geometricRegion(const dictionary& dict);
+
+	~geometricRegion() = default;
+
+	//// - methods
+	bool   isInside(const realx3& p) const;
+
+	realx3 peek() const;
+
+	//// IO operation
+	bool   read(const dictionary& dict);
+
+	bool   write(dictionary& dict) const;
+};
+
+} // namespace pFlow
+
+#include "geometricRegion.cpp"
+
+#endif

src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegions.cpp

@@ -0,0 +1,9 @@
+
+
+
+#include "geometricRegions.hpp"
+
+
+template class pFlow::geometricRegion<pFlow::box>;
+
+template class pFlow::geometricRegion<pFlow::sphere>;

src/phasicFlow/structuredData/peakableRegion/geometricRegion/geometricRegions.hpp

@@ -18,30 +18,23 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
+#ifndef __geometricRegions_hpp__
+#define __geometricRegions_hpp__
 
-#ifndef __demParticles_hpp__ 
+#include "geometricRegion.hpp"
-#define __demParticles_hpp__
+#include "box.hpp"
-
+#include "sphere.hpp"
-#include "demComponent.hpp"
+#include "cylinder.hpp"
 
 namespace pFlow
 {
 
+using boxRegion = geometricRegion<box>;
 
-class demParticles
+using sphereRegion = geometricRegion<sphere>;
-:
-	public demComponent
-{
 
-public:
+using cylinderRegion = geometricRegion<cylinder>;
-
-	demParticles(systemControl& control):
-		demComponent("particles", control)
-	{}
-
-
-};
 
 }
 
-#endif
+#endif //__geometricRegions_hpp__

src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.cpp

@@ -18,13 +18,13 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #include "peakableRegion.hpp"
 #include "dictionary.hpp"
 
 pFlow::peakableRegion::peakableRegion
 (
-	const word& type, const dictionary& dict
+	const word& type, 
+	const dictionary& dict
 )
 {
 	CONSUME_PARAM(type);
@@ -38,7 +38,7 @@ pFlow::uniquePtr<pFlow::peakableRegion> pFlow::peakableRegion::create
 	const dictionary& dict
 )
 { 
-	word regionType = angleBracketsNames("peakableRegion", type);
+	word regionType = angleBracketsNames("peakable", type);
 
 	if( wordvCtorSelector_.search(regionType) )
 	{
@@ -48,7 +48,7 @@ pFlow::uniquePtr<pFlow::peakableRegion> pFlow::peakableRegion::create
 	{
 		printKeys
 		( 
-			fatalError << "Ctor Selector "<< regionType << " dose not exist. \n"
+			fatalError << "Ctor Selector "<< regionType << " does not exist. \n"
 			<<"Avaiable ones are: \n\n"
 			,
 			wordvCtorSelector_

src/phasicFlow/structuredData/peakableRegion/peakableRegion/peakableRegion.hpp

@@ -18,25 +18,19 @@ Licence:
 
 -----------------------------------------------------------------------------*/
 
-
 #ifndef __peakableRegion_hpp__
 #define __peakableRegion_hpp__
 
 #include "types.hpp"
 #include "virtualConstructor.hpp"
 
-
 namespace pFlow
 {
 
-
 class dictionary;
 
 class peakableRegion
 {
-protected:
-
-	
 public:
 
 	// type info
@@ -63,7 +57,6 @@ public:
 
 	virtual realx3 peek() const = 0;
 
-
 	//// - IO operatoins
 
 	virtual bool   read(const dictionary& dict) = 0;
@@ -73,10 +66,8 @@ public:
 	// - static create
 	static uniquePtr<peakableRegion>
 	create(const word& type, const dictionary& dict);
-
 };
 
-} // pFlow
+} // namespace pFlow
-
 
 #endif

src/phasicFlow/structuredData/peakableRegion/peakableRegionInstantiate.cpp

@@ -1,30 +0,0 @@
-/*------------------------------- phasicFlow ---------------------------------
-      O        C enter of
-     O O       E ngineering and
-    O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
-------------------------------------------------------------------------------
-  Copyright (C): www.cemf.ir
-  email: hamid.r.norouzi AT gmail.com
-------------------------------------------------------------------------------  
-Licence:
-  This file is part of phasicFlow code. It is a free software for simulating 
-  granular and multiphase flows. You can redistribute it and/or modify it under
-  the terms of GNU General Public License v3 or any other later versions. 
- 
-  phasicFlow is distributed to help others in their research in the field of 
-  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
-  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-
------------------------------------------------------------------------------*/
-
-#include "boxRegion.hpp"
-#include "sphereRegion.hpp"
-#include "cylinderRegion.hpp"
-
-
-template class pFlow::PeakableRegion<pFlow::boxRegion>;
-
-template class pFlow::PeakableRegion<pFlow::sphereRegion>;
-
-template class pFlow::PeakableRegion<pFlow::cylinderRegion>;

src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp

@@ -193,14 +193,14 @@ const pFlow::pFlagTypeHost&
 	return pFlagsH_;
 }
 
-FUNCTION_H
+
 const typename pFlow::internalPoints::PointsType& 
     pFlow::internalPoints::pointPosition()const
 {
     return pointPosition_;
 }
 
-FUNCTION_H
+
 typename pFlow::internalPoints::PointsType&
 	pFlow::internalPoints::pointPosition()
 {
@@ -255,7 +255,7 @@ bool pFlow::internalPoints::deletePoints
     return true;
 }
 
-FUNCTION_H
+
 pFlow::uint32 pFlow::internalPoints::updateFlag
 (
 	const domain& dm, 
@@ -296,7 +296,108 @@ void pFlow::internalPoints::fillNeighborsLists
 		
 }
 
-FUNCTION_H
+bool
+pFlow::internalPoints::insertPoints(
+  const realx3Vector& points,
+  anyList&            varList
+)
+{
+	uint32 numNew = static_cast<uint32>(points.size());
+
+	auto aRange = pFlagsD_.activeRange();
+	uint32 emptySpots = pFlagsD_.capacity() - pFlagsD_.numActive();
+	message msg;
+
+	if( numNew > emptySpots )
+	{
+		// increase the capacity to hold new points 
+		aRange = pFlagsD_.activeRange();
+		uint32 newCap = pFlagsD_.changeCapacity(numNew);
+		unSyncFlag();
+		varList.emplaceBack(
+			msg.addAndName(message::CAP_CHANGED),
+			newCap);
+	}
+	
+	
+	// first check if it is possible to add to the beggining of list
+	if(numNew <= aRange.start())
+	{
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			0u, numNew);
+	}
+	// check if it is possible to add to the end of the list  
+	else if( numNew <= pFlagsD_.capacity() - aRange.end() )
+	{	 	
+		varList.emplaceBack<uint32IndexContainer>(
+			msg.addAndName(message::ITEM_INSERT),
+			aRange.end(), aRange.end()+numNew);
+	}
+	// we should fill the scattered empty spots
+	else
+	{
+		notImplementedFunction;
+		return false;
+	}
+
+	const auto& indices = varList.getObject<uint32IndexContainer>(
+		message::eventName(message::ITEM_INSERT)
+	);
+
+	auto nAdded = pFlagsD_.addInternalPoints(indices.deviceView());
+	unSyncFlag();
+
+	if(nAdded != numNew )
+	{
+		fatalErrorInFunction;
+		return false;
+	}
+
+	pointPosition_.reserve( pFlagsD_.capacity() );
+	if(!pointPosition_.insertSetElement(indices, points))
+	{
+		fatalErrorInFunction<<
+		"Error in inserting new positions into pointPosition"<<
+		" internalPoints field"<<endl;
+		return false;
+	}
+
+	auto newARange = pFlagsD_.activeRange();
+
+	if( aRange.end() != newARange.end() )
+	{
+		varList.emplaceBack(
+			msg.addAndName(message::RANGE_CHANGED),
+			newARange);
+		varList.emplaceBack(
+			msg.addAndName(message::SIZE_CHANGED),
+			newARange.end());
+	}
+	else if(aRange.start() != newARange.start())
+	{
+		varList.emplaceBack(
+			msg.addAndName(message::RANGE_CHANGED),
+			newARange);
+	}
+	
+	auto tInfo = time().TimeInfo();
+
+	if(!notify(
+		tInfo.iter(),
+		tInfo.t(),
+		tInfo.dt(),
+		msg, 
+		varList))
+	{
+		fatalErrorInFunction<<
+		"Error in notifying observers of "<< subscriberName()<<endl;
+		return false;
+	}
+
+	return true;
+}
+
 bool pFlow::internalPoints::read
 (
 	iIstream& is
@@ -322,7 +423,7 @@ bool pFlow::internalPoints::read
 }
 
 
-FUNCTION_H
+
 bool pFlow::internalPoints::write
 (
 	iOstream& os
@@ -332,7 +433,7 @@ bool pFlow::internalPoints::write
 	return aPoints.write(os);
 }
 
-FUNCTION_H
+
 bool pFlow::internalPoints::read
 (
 	iIstream& is, 
@@ -359,7 +460,7 @@ bool pFlow::internalPoints::read
 		
 
 
-FUNCTION_H
+
 bool pFlow::internalPoints::write
 (
 	iOstream& os, 

src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.hpp

@@ -83,6 +83,7 @@ private:
 			uint32 fromBoundaryIndex,
 			uint32 toBoundaryIndex);
 
+
 protected:
 	
 	void syncPFlag()const;
@@ -134,27 +135,25 @@ public:
 
 
 	//// - Methods
-		FUNCTION_H
+		
 		const pFlagTypeDevice& activePointsMaskDevice()const;
 
 		
-		FUNCTION_H
 		const pFlagTypeHost& activePointsMaskHost()const;
 
 		// - Const access pointPosition
-		FUNCTION_H
 		const PointsType& pointPosition()const;
 		
-		FUNCTION_H
+		
         PointsType& pointPosition();
 
-        INLINE_FUNCTION_H
+        inline
         auto pointPositionHost()const
         {
             return pointPosition_.hostView();
         }
 
-        INLINE_FUNCTION_H
+        inline
         auto pointPositionDevice()const
         {
             return pointPosition_.deviceView();
@@ -163,40 +162,40 @@ public:
 		PointsTypeHost activePointsHost()const;
 
 		// - size of data structure
-		INLINE_FUNCTION_H
+		inline
 		uint32 size()const
 		{
 			return pFlagsD_.activeRange_.end();
 		}
 		
 		// - maximum capacity of data structure 
-		INLINE_FUNCTION_H
+		inline
 		uint32 capacity()const
 		{
-			return pointPosition_.capacity();
+			return pFlagsD_.capacity();
 		}
 
-		INLINE_FUNCTION_H
+		inline
 		bool empty()const
 		{
 			return size()==0u;
 		}
 		
 		// - number of active points 
-		INLINE_FUNCTION_H
+		inline
 		uint32 numActive() const
 		{
 			return pFlagsD_.numActive();
 		}
 				
 		// - if all points are active 
-		INLINE_FUNCTION_H
+		inline
 		bool isAllActive()const
 		{
 			return pFlagsD_.isAllActive();
 		}
 
-		INLINE_FUNCTION_H
+		inline
 		auto activeRange()const
 		{
 			return pFlagsD_.activeRange();
@@ -216,16 +215,16 @@ public:
 
 		///@brief delete points at indices given in delPoints.
 		/// The default is that delPoints contains sorted indices 
-		FUNCTION_H
+		
 		bool deletePoints(
 			scatteredFieldAccess<uint32, memory_space> delPoints);
 		
-		FUNCTION_H
+		
 		uint32 updateFlag(
 			const domain& dm, 
 			const std::array<real,6>& dist);
 		
-		FUNCTION_H
+		
 		void fillNeighborsLists(
 			ViewType1D<uint32, memory_space> leftList,
 			ViewType1D<uint32, memory_space> rightList,
@@ -235,26 +234,23 @@ public:
 			ViewType1D<uint32, memory_space> frontList);
 		
 		
+		bool insertPoints(const realx3Vector& points, anyList& varList);
 		
 	
 	//// - IO operations 
 		
 		/// Read 
-		FUNCTION_H
 		bool read(iIstream& is);
 		
 
 		/// Write 
-		FUNCTION_H
 		bool write(iOstream& os)const;
 
 		/// Read 
-		FUNCTION_H
 		bool read(iIstream& is, const IOPattern& iop);
 		
 
 		/// Write 
-		FUNCTION_H
 		bool write(iOstream& os, const IOPattern& iop)const;	
 
 };

src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlag.hpp

@@ -20,10 +20,10 @@ Licence:
 #ifndef __pointFlag_hpp__
 #define __pointFlag_hpp__
 
-
 #include "phasicFlowKokkos.hpp"
 #include "domain.hpp"
 #include "scatteredFieldAccess.hpp"
+#include "globalSettings.hpp"
 
 namespace pFlow
 {
@@ -112,8 +112,7 @@ public:
 
 	friend class internalPoints;
 
-	pointFlag()
+	pointFlag() = default;
-	{}
 
 	pointFlag(uint32 capacity, uint32 start, uint32 end )
 	:
@@ -122,7 +121,7 @@ public:
 		activeRange_(start, end),
 		isAllActive_(true)
 	{
-		fill(flags_, 0, capacity, static_cast<uint8>(Flag::DELETED));
+		fill(flags_, end, capacity, static_cast<uint8>(Flag::DELETED));
 		fill(flags_, activeRange_, static_cast<uint8>(Flag::INTERNAL));
 	}
 
@@ -160,6 +159,12 @@ public:
 		return activeRange_;
 	}
 
+	INLINE_FUNCTION_HD
+	uint32 capacity()const
+	{
+		return flags_.extent(0);
+	}
+
 	INLINE_FUNCTION_HD
 	auto numActive()const
 	{
@@ -289,6 +294,9 @@ public:
 		const box& validBox,
 		ViewType1D<realx3, memory_space> points);
 
+	uint32 addInternalPoints(
+		const ViewType1D<uint32, memory_space>& points);
+
 	bool deletePoints(
 		scatteredFieldAccess<uint32, memory_space> points);
 
@@ -345,6 +353,8 @@ public:
 		ViewType1D<uint32, memory_space> 	indices,
 		ViewType1D<uint32, memory_space> 	onOff);
 
+	uint32 changeCapacity( uint32 reqEmptySpots );
+
 
 
 };

src/phasicFlow/structuredData/pointStructure/internalPoints/pointFlagKernels.hpp

@@ -122,62 +122,52 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markOutOfBoxDelete
 	return numDeleted;
 }
 
-/*template<typename ExecutionSpace>
+template<typename ExecutionSpace>
-pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::scanPointFlag()
+pFlow::uint32 
+pFlow::pointFlag<ExecutionSpace>::addInternalPoints
+(
+	const ViewType1D<uint32, memory_space>& points
+)
 {
+	if(points.extent(0)==0u)return 0;
 
-	using rpScanFlag = Kokkos::RangePolicy<execution_space,  
+	uint32 minRange;
-			Kokkos::IndexType<uint32>>;
+	uint32 maxRange;
-
+	uint32 numAdded = 0;
-	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",
+		"pointFlagKernels::addInternalPoints",
-			rpScanFlag(start, end),
+		deviceRPolicyStatic(0, points.extent(0)),
-			CLASS_LAMBDA_HD(
+		LAMBDA_HD(
 			uint32 i, 
 			uint32& minUpdate, 
 			uint32& maxUpdate, 
-				uint32& sumToUpdate)
+			uint32& addToUpdate)
 		{
-				if(this->isActive(i))
+			uint32 idx = points(i);
-				{
+			if( flags_[idx] <= DELETED) addToUpdate ++;
-					sumToUpdate++;
+			minUpdate = min(minUpdate,idx);
-					minUpdate = min(minUpdate,i);
+			maxUpdate = max(maxUpdate,idx);
-					maxUpdate = max(maxUpdate,i);
+			flags_[idx] = INTERNAL;	
-				}
 		},
 		Kokkos::Min<uint32>(minRange), 
 		Kokkos::Max<uint32>(maxRange),
-			numActive);	
+		numAdded);	
-	}
 
-	if(numActive==0)
+	// add one to max range to make range it half open 
-	{
-		minRange = 0;
-		maxRange = 0;
-	}
-	else
-	{
-		// add one to maxRange to make it half-open
 	maxRange++;
-	}
+	
+	minRange = min(activeRange_.start(), minRange);
+	maxRange = max(activeRange_.end(), maxRange);
 
 	activeRange_ = {minRange, maxRange};
-	numActive_ 	= numActive;
+	numActive_ 	= numActive_ + numAdded;
+	
 	isAllActive_ = activeRange_.numElements() == numActive_;
 	
-	return numActive;
+	return numAdded;
-}*/
+}
 
 template<typename ExecutionSpace>
 bool pFlow::pointFlag<ExecutionSpace>::deletePoints
@@ -187,12 +177,6 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 {
 	if(points.empty())return true;
 	
-	//uint32 minIndex = points.getFirstCopy();
-	//uint32 maxIndex = points.getLastCopy();
-	
-	//if(maxIndex<minIndex) return false;
-	//if(maxIndex>activeRange_.end())return false;
-	//if(minIndex<activeRange_.start())return false;
 	
 	using policy = Kokkos::RangePolicy<
 		execution_space,  
@@ -237,10 +221,6 @@ bool pFlow::pointFlag<ExecutionSpace>::deletePoints
 	uint32 s = points.size();
 	if(s==0u)return true;
 	
-	//if(maxIndex<minIndex) return false;
-	//if(maxIndex>activeRange_.end())return false;
-	//if(minIndex<activeRange_.start())return false;
-
 	using policy = Kokkos::RangePolicy<
 		execution_space,  
 		Kokkos::IndexType<uint32>>;
@@ -549,4 +529,92 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markDelete
 	return numMarked;
 }
 
+template<typename ExecutionSpace>
+pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::changeCapacity
+(
+	uint32 reqEmptySpots
+)
+{
+	
+	uint32 oldCap = capacity();
+	uint32 newCap = oldCap* pFlow::gSettings::vectorGrowthFactor__+1;
+	uint32 emptySpots = newCap - numActive_;
+
+	while( emptySpots < reqEmptySpots )
+	{
+		newCap = newCap*pFlow::gSettings::vectorGrowthFactor__+1;
+		uint32 emptySpots = newCap - numActive_;
+	}
+	
+	viewType newFlags(flags_.label(), newCap);
+	// copy contnent 
+	copy(newFlags, 0u, flags_, 0u, oldCap);
+	fill(newFlags, oldCap, newCap, static_cast<uint8>(DELETED));
+
+	// reference replacement 
+	flags_ = newFlags;
+
+	return newCap;
+}
+
 #endif // __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;
+}*/