boundaryExit beforeIteration, not tested

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -274,7 +274,7 @@ pFlow::sphereParticles::sphereParticles(
 		(
 			"propertyId",
 			"",
-			objectFile::READ_ALWAYS,
+			objectFile::READ_NEVER,
 			objectFile::WRITE_NEVER
 		),
 		dynPointStruct(),

src/Particles/particles/particles.hpp

@@ -87,11 +87,11 @@ protected:
 		return dynPointStruct_.pointPosition();	
 	}
 
-	inline 
+	/*inline 
 	auto& velocity()
 	{
 		return dynPointStruct_.velocity();
-	}
+	}*/
 
 	inline 
 	auto& shapeIndex()

src/phasicFlow/CMakeLists.txt

@@ -76,7 +76,7 @@ structuredData/pointStructure/pointStructure.cpp
 structuredData/boundaries/boundaryBase/boundaryBase.cpp
 structuredData/boundaries/boundaryExit/boundaryExit.cpp
 structuredData/boundaries/boundaryNone/boundaryNone.cpp
-structuredData/pointStructure/boundaryList.cpp
+structuredData/boundaries/boundaryList.cpp
 structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp
 structuredData/pointStructure/selectors/selectBox/selectBox.cpp
 structuredData/pointStructure/selectors/selectRange/selectRange.cpp 

src/phasicFlow/Kokkos/ViewAlgorithms.hpp

@@ -241,19 +241,19 @@ void copy(
 }
 
 template <
-	typename dType,
+	typename Type,
-	typename sType,
 	typename... sProperties>
 INLINE_FUNCTION_H
 void getNth(
-	dType& dst,
+	Type& dst,
-	const ViewType1D<sType, sProperties...>& src,
+	const ViewType1D<Type, sProperties...>& src,
 	const uint32 n
 	)
 {
-	range32 span(n,n+1);
+		
-	auto subV = Kokkos::subview(src, span);
+	auto subV = Kokkos::subview(src, Kokkos::make_pair(n,n+1));
-	hostViewType1D<dType> dstView("getNth",1);
+	hostViewType1D<Type> dstView("getNth",1);
+	//hostViewTypeScalar
 	Kokkos::deep_copy(dstView,subV);
 	dst = *dstView.data();
 }
@@ -425,13 +425,12 @@ int32 binarySearch(
 template<
 	typename Type,
 	typename... properties,
-	typename dType,
 	typename... dProperties>
 void exclusiveScan(
 	const ViewType1D<Type, properties...>& view,
 	uint32 start,
 	uint32 end,
-	ViewType1D<dType, dProperties...>& dView,
+	ViewType1D<Type, dProperties...>& dView,
 	uint32 dStart )
 {
 
@@ -439,7 +438,7 @@ void exclusiveScan(
 	(
 		areAccessible<
 			typename ViewType1D<Type, properties...>::execution_space,
-			typename ViewType1D<dType, dProperties...>::memory_space>(),
+			typename ViewType1D<Type, dProperties...>::memory_space>(),
 		"In exclusiveScan, view and dView should have the same space"
 
 	);
@@ -448,7 +447,7 @@ void exclusiveScan(
 	
 	uint32 numElems = end-start;
 	
-	pFlow::algorithms::KOKKOS::exclusiveScan<Type,dType,ExecutionSpace>(
+	pFlow::algorithms::KOKKOS::exclusiveScan<Type,ExecutionSpace>(
 		view.data()+start,
 		dView.data()+dStart,
 		numElems);
@@ -458,13 +457,12 @@ void exclusiveScan(
 template<
 	typename Type,
 	typename... properties,
-	typename dType,
 	typename... dProperties>
 void inclusiveScan(
 	const ViewType1D<Type, properties...>& view,
 	uint32 start,
 	uint32 end,
-	ViewType1D<dType, dProperties...>& dView,
+	ViewType1D<Type, dProperties...>& dView,
 	uint32 dStart)
 {
 	using ExecutionSpace = typename ViewType1D<Type, properties...>::execution_space;
@@ -473,14 +471,14 @@ void inclusiveScan(
 	(
 		areAccessible<
 			typename ViewType1D<Type, properties...>::execution_space,
-			typename ViewType1D<dType, dProperties...>::memory_space>(),
+			typename ViewType1D<Type, dProperties...>::memory_space>(),
 		"In exclusiveScan, view and dView should have the same space"
 	);
 
 
 	uint32 numElems = end-start;
 
-	pFlow::algorithms::KOKKOS::inclusiveScan<Type,dType,ExecutionSpace>(
+	pFlow::algorithms::KOKKOS::inclusiveScan<Type,ExecutionSpace>(
 		view.data()+start,
 		dView.data()+dStart,
 		numElems);

src/phasicFlow/algorithms/kokkosAlgorithms.hpp

@@ -30,15 +30,15 @@ namespace pFlow::algorithms::KOKKOS
 
 template<typename Type, typename ExecutionSpace>
 INLINE_FUNCTION_H
-int32 count(const Type* first, int32 numElems, const Type& val)
+uint32 count(const Type* first, uint32 numElems, const Type& val)
 {
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
-	int32 num = 0;
+	uint32 num = 0;
 	Kokkos::parallel_reduce("count",
 							policy(0, numElems),
-							LAMBDA_HD(int32 i, int32& updateVal){
+							LAMBDA_HD(uint32 i, uint32& updateVal){
 								if(equal(first[i],val)) updateVal++;
 							},
 							num);
@@ -50,17 +50,17 @@ int32 count(const Type* first, int32 numElems, const Type& val)
 
 template<typename Type, typename ExecutionSpace>
 INLINE_FUNCTION_H
-void fillSequence(Type* first, int32 numElems, const Type& firstVal)
+void fillSequence(Type* first, uint32 numElems, const Type& firstVal)
 {
 	
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 
 	Kokkos::parallel_for(
 		"fillSequence",
 		policy(0, numElems),
-		LAMBDA_HD(int32 i){ 
+		LAMBDA_HD(uint32 i){ 
 			first[i] = firstVal+i;
 		});
 	Kokkos::fence();
@@ -68,15 +68,15 @@ void fillSequence(Type* first, int32 numElems, const Type& firstVal)
 
 template<typename Type, typename indexType, typename ExecutionSpace>
 INLINE_FUNCTION_H
-void fillSelected(Type* first, const indexType* indices, const int32 numElems, const Type val)
+void fillSelected(Type* first, const indexType* indices, const uint32 numElems, const Type val)
 {
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 	Kokkos::parallel_for(
 		"fillSelected",
 		policy(0,numElems),
-		LAMBDA_HD(int32 i){
+		LAMBDA_HD(uint32 i){
 			first[indices[i]]= val;
 		});
 	Kokkos::fence();
@@ -84,16 +84,16 @@ void fillSelected(Type* first, const indexType* indices, const int32 numElems, c
 
 template<typename Type, typename indexType, typename ExecutionSpace>
 INLINE_FUNCTION_H
-void fillSelected(Type* first, const indexType* indices, const Type* vals, const int32 numElems)
+void fillSelected(Type* first, const indexType* indices, const Type* vals, const uint32 numElems)
 {
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 
 	Kokkos::parallel_for(
 		"fillSelected",
 		policy(0,numElems),
-		LAMBDA_HD(int32 i){
+		LAMBDA_HD(uint32 i){
 			first[indices[i]]= vals[i];
 		});
 	Kokkos::fence();
@@ -101,17 +101,17 @@ void fillSelected(Type* first, const indexType* indices, const Type* vals, const
 
 template<typename Type, typename ExecutionSpace>
 INLINE_FUNCTION_H
-Type max(const Type* first, int32 numElems)
+Type max(const Type* first, uint32 numElems)
 {
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 	Type maxElement=0;
 
 	Kokkos::parallel_reduce(
 		"max",
 		policy(0, numElems),
-		LAMBDA_HD(int32 i, Type& maxUpdate){
+		LAMBDA_HD(uint32 i, Type& maxUpdate){
 			if(maxUpdate<first[i]) maxUpdate = first[i];
 		},
 		Kokkos::Max<Type>(maxElement));
@@ -144,38 +144,38 @@ Type min(const Type* first, int32 numElems)
 //void sort(Type* first, int32 numElems, CompareFunc compare);
 //void permuteSort(const Type* first, PermuteType* pFirst, int32 numElems);
 
-template<typename Type, typename DestType, typename ExecutionSpace>
+template<typename Type, typename ExecutionSpace>
-void exclusiveScan(Type* first, DestType* dFirst, int32 numElems)
+void exclusiveScan(Type* first, Type* dFirst, uint32 numElems)
 {	
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 
 	Kokkos::parallel_scan(
 		"exclusiveScan",
 		policy(0, numElems),
-		LAMBDA_HD(const int32 i, DestType& valToUpdate, const bool final)
+		LAMBDA_HD(const uint32 i, Type& valToUpdate, const bool final)
 		{
-			const int32 val = first[i];
+			const Type val = first[i];
 			if(final)
 				dFirst[i] = valToUpdate;
 			valToUpdate += val;
 		});
 }
 
-template<typename Type, typename DestType, typename ExecutionSpace>
+template<typename Type, typename ExecutionSpace>
-void inclusiveScan(Type* first, DestType* dFirst, int32 numElems)
+void inclusiveScan(Type* first, Type* dFirst, uint32 numElems)
 {
 	using policy = Kokkos::RangePolicy<
 			ExecutionSpace,
-			Kokkos::IndexType<int32> >;
+			Kokkos::IndexType<uint32> >;
 
 	Kokkos::parallel_scan(
 		"inclusiveScan",
 		policy(0, numElems),
-		LAMBDA_HD(const int32 i, int32& valToUpdate, const bool final)
+		LAMBDA_HD(const uint32 i, Type& valToUpdate, const bool final)
 		{
-			const int32 val = first[i];
+			const Type val = first[i];
 			valToUpdate += val;
 			if(final)
 				dFirst[i] = valToUpdate;

src/phasicFlow/containers/Vector/stdVectorHelper.hpp

@@ -48,6 +48,23 @@ public:
 template<typename T>
 using vecAllocator = std::allocator<T>;
 
+template<typename T>
+inline
+span<T> makeSpan(std::vector<T>& container)
+{
+    return span<T>(container.data(), container.size());
+}
+
+template<typename T>
+inline
+span<T> makeSpan(const std::vector<T>& container)
+{
+    return span<T>(
+		const_cast<T*>(container.data()), 
+		container.size());
+}
+
+
 template<typename T>
 inline
 bool writeSpan(

src/phasicFlow/containers/pointField/pointField.cpp

@@ -1,4 +1,3 @@
-#include "pointField.hpp"
 /*------------------------------- phasicFlow ---------------------------------
       O        C enter of
      O O       E ngineering and
@@ -41,7 +40,25 @@ pFlow::pointField<T, MemorySpace>::pointField
 	boundaryFieldList_(pStruct.boundaries(), *this),
 	pStruct_(pStruct),
 	defaultValue_(defVal)
-{}
+{
+    
+    if(IOobject::implyRead())
+    {
+        REPORT(1)<< "Reading field "<< Green_Text(IOobject::name())<<
+        " from "<<IOobject::path()<<END_REPORT;
+    }
+    else
+    {
+        REPORT(1)<< "Creating field "<< Green_Text(IOobject::name())<<END_REPORT;
+    }
+    if( !IOobject::readObject() )
+    {
+        fatalErrorInFunction<<
+        "Error in reading from file "<<IOobject::path()<<endl;
+        fatalExit;
+    }    
+
+}
 
 template<class T, class MemorySpace>
 pFlow::pointField<T, MemorySpace>::pointField
@@ -66,7 +83,23 @@ pFlow::pointField<T, MemorySpace>::pointField
 	boundaryFieldList_(pStruct.boundaries(), *this),
 	pStruct_(pStruct),
 	defaultValue_(defVal)
-{}
+{
+    if(IOobject::implyRead())
+    {
+        REPORT(1)<< "Reading field "<< Green_Text(IOobject::name())<<
+        " from "<<IOobject::path()<<END_REPORT;
+    }
+    else
+    {
+        REPORT(1)<< "Creating field "<< Green_Text(IOobject::name())<<END_REPORT;
+    }
+    if( !IOobject::readObject() )
+    {
+        fatalErrorInFunction<<
+        "Error in reading from file "<<IOobject::path()<<endl;
+        fatalExit;
+    }
+}
 
 template<class T, class MemorySpace>
 pFlow::pointField<T, MemorySpace>::pointField
@@ -98,7 +131,9 @@ pFlow::pointField<T, MemorySpace>::pointField
 	boundaryFieldList_(pStruct.boundaries(), *this),
 	pStruct_(pStruct),
 	defaultValue_(defVal)
-{}
+{
+
+}
 
 
  /*

src/phasicFlow/containers/pointField/pointField.hpp

@@ -52,7 +52,7 @@ public:
 
 	using execution_space 	= typename InternalFieldType::execution_space;
 
-protected:
+private:
 
     //// - data members 
 

src/phasicFlow/containers/span/span.hpp

@@ -157,7 +157,7 @@ public:
 
 };
 
-template<typename T, typename... properties, template<class, class...> class Container>
+/*template<typename T, typename... properties, template<class, class...> class Container>
 size_t makeSpan(Container<T, properties...>& container)
 {
     return span<T>(container.data(), container.size());
@@ -169,7 +169,7 @@ size_t makeSpan(const Container<T, properties...>& container)
     return span<T>(
         const_cast<T*>(container.data()), 
         container.size());
-}
+}*/
 
 
 template<typename T>
@@ -192,7 +192,7 @@ span<const char> charSpan(span<const T> s)
         s.size()*el);
 }
 
-template<typename T, template<class> class Container>
+/*template<typename T, template<class> class Container>
 span<T> makeSpan(Container<T>& container)
 {
     return span<T>(container.data(), container.size());
@@ -204,7 +204,7 @@ span<T> makeSpan(const Container<T>& container)
     return span<T>(
         const_cast<T*>(container.data()), 
         container.size());
-}
+}*/
 
 
 

src/phasicFlow/repository/Time/baseTimeControl.cpp

@@ -25,7 +25,8 @@ Licence:
 pFlow::baseTimeControl::baseTimeControl
 (
     const dictionary &dict,
-    const word& intervalPrefix
+    const word& intervalPrefix,
+	real defStartTime
 )
 {
 	auto tControl = dict.getVal<word>("timeControl");
@@ -47,7 +48,7 @@ pFlow::baseTimeControl::baseTimeControl
 
 	if(isTimeStep_)
 	{
-		auto startTime = (dict.getValOrSet<real>("startTime", 0.0));
+		auto startTime = (dict.getValOrSet<real>("startTime", defStartTime));
     	auto endTime = (dict.getValOrSet<real>("endTime", largeValue));
 		auto interval = dict.getVal<real>(intervalWord);
 		rRange_ = realStridedRange(startTime, endTime, interval);
@@ -55,7 +56,7 @@ pFlow::baseTimeControl::baseTimeControl
 	}
 	else
 	{
-		auto startTime = (dict.getValOrSet<int32>("startTime", 0.0));
+		auto startTime = (dict.getValOrSet<int32>("startTime", 0));
     	auto endTime = (dict.getValOrSet<int32>("endTime", largestPosInt32));
 		auto interval = dict.getVal<int32>(intervalWord);
 		iRange_ = int32StridedRagne(startTime, endTime, interval);

src/phasicFlow/repository/Time/baseTimeControl.hpp

@@ -1,20 +1,20 @@
 /*------------------------------- phasicFlow ---------------------------------
-      O        C enter of
+			O        C enter of
-     O O       E ngineering and
+		 O O       E ngineering and
-    O   O      M ultiscale modeling of
+		O   O      M ultiscale modeling of
-   OOOOOOO     F luid flow       
+	 OOOOOOO     F luid flow       
 ------------------------------------------------------------------------------
-  Copyright (C): www.cemf.ir
+	Copyright (C): www.cemf.ir
-  email: hamid.r.norouzi AT gmail.com
+	email: hamid.r.norouzi AT gmail.com
 ------------------------------------------------------------------------------  
 Licence:
-  This file is part of phasicFlow code. It is a free software for simulating 
+	This file is part of phasicFlow code. It is a free software for simulating 
-  granular and multiphase flows. You can redistribute it and/or modify it under
+	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. 
+	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 field of 
-  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+	granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
-  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+	implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 -----------------------------------------------------------------------------*/
 #ifndef __baseTimeControl_hpp__
@@ -30,8 +30,8 @@ namespace pFlow
 class baseTimeControl
 {
 private:
-   
+	 
-    bool 	isTimeStep_;
+	bool 	isTimeStep_;
 
 	int32StridedRagne iRange_;
 
@@ -40,7 +40,10 @@ private:
 
 public:
 
-	baseTimeControl(const dictionary& dict, const word& intervalPrefix = "");
+	baseTimeControl(
+		const dictionary& dict, 
+		const word& intervalPrefix = "",
+		real defStartTime = 0.0);
 
 
 	bool timeEvent(uint32 iter, real t, real dt)const;

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

@@ -50,9 +50,9 @@ pFlow::boundaryBase::boundaryBase
 :
 	subscriber(dict.name()),
 	boundaryPlane_(bplane),
+	indexList_(groupNames(dict.name(),"indexList")),
 	neighborLength_(dict.getVal<real>("neighborLength")),
 	internal_(internal),
-	indexList_(groupNames(dict.name(),"indexList")),
 	mirrorProcessoNo_(dict.getVal<uint32>("mirrorProcessorNo")),
 	name_(dict.name()),
 	type_(dict.getVal<word>("type"))

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

@@ -25,7 +25,7 @@ Licence:
 #include "subscriber.hpp"
 #include "VectorSingles.hpp"
 #include "plane.hpp"
-#include "scatterFieldAccess.hpp"
+#include "scatteredFieldAccess.hpp"
 
 #include "streams.hpp"
 
@@ -44,29 +44,27 @@ class boundaryBase
 public:
 
     using pointFieldAccessType = 
-        scatterFieldAccess<realx3,DefaultExecutionSpace>;
+        deviceScatteredFieldAccess<realx3>;
 
-
+private:
-protected:
 
 	const plane& 	boundaryPlane_;
 
+	/// list of particles indices on device 
+	uint32Vector_D 	indexList_;
+
 	/// The length defined for creating neighbor list 
 	real  			neighborLength_;	
 
 	/// a reference to 
 	internalPoints& internal_;
 
-	/// list of particles indices on device 
-	uint32Vector_D 	indexList_;
-
 	uint32 			mirrorProcessoNo_;
 
 	word 		name_;
 
 	word 		type_;
 
-
 public:
 
 	TypeInfo("boundaryBase");
@@ -85,7 +83,7 @@ public:
 
 	boundaryBase& operator=(boundaryBase&&) = default;
 
-	virtual ~boundaryBase() = default;
+	~boundaryBase() override = default;
 	
 
 	create_vCtor
@@ -117,22 +115,42 @@ public:
 		return name_;
 	}
 
+	bool empty()const
+	{
+		return indexList_.size()==0;
+	}
+
 	auto size()const
 	{
 		return indexList_.size();
 	}
 
+	const plane& boundaryPlane()const
+	{
+		return boundaryPlane_;
+	}
+
 	auto capacity()const
 	{
 		return indexList_.capacity();
 	}
 
+	const internalPoints& internal()const
+	{
+		return internal_;
+	}
+
+	internalPoints& internal()
+	{
+		return internal_;
+	}
+
 	/// @brief set the size of indexList 
 	/// Always make sure that size+1 <= capacity
 	void setSize(uint32 newSize);
 
 	virtual 
-    bool beforeIteratoin(uint32 iterNum, real t) = 0 ;
+    bool beforeIteration(uint32 iterNum, real t, real dt) = 0 ;
 
 	virtual 
     bool iterate(uint32 iterNum, real t) = 0;
@@ -146,6 +164,11 @@ public:
 		return indexList_;
 	}
 
+	auto& indexList()
+	{
+		return indexList_;
+	}
+
     pointFieldAccessType thisPoints();
 
     virtual

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

@@ -1,6 +1,6 @@
 
-#ifndef __scatterFieldAccess_hpp__
+#ifndef __scatteredFieldAccess_hpp__
-#define __scatterFieldAccess_hpp__
+#define __scatteredFieldAccess_hpp__
 
 
 #include "phasicFlowKokkos.hpp"
@@ -9,8 +9,8 @@ namespace pFlow
 {
 
 
-template<typename T, typename MemorySpace>
+template<typename T, typename MemorySpace=void>
-class scatterFieldAccess
+class scatteredFieldAccess
 {
 public:
 	
@@ -22,7 +22,7 @@ public:
 
     using execution_space 	= typename viewType::execution_space;
 
-protected:
+private:
 
 	uint32 								size_ = 0;
 
@@ -32,12 +32,12 @@ protected:
 
 public:
 
-    scatterFieldAccess():
+    scatteredFieldAccess():
         indices_(),
         fieldVals_()    
     {}
 
-	scatterFieldAccess(
+	scatteredFieldAccess(
 		uint32 	sz,
 		ViewType1D<uint32, memory_space> ind,
 		ViewType1D<T, memory_space> fVals)
@@ -47,15 +47,15 @@ public:
 		fieldVals_(fVals)
 	{}
 
-	scatterFieldAccess(const scatterFieldAccess&) = default;
+	scatteredFieldAccess(const scatteredFieldAccess&) = default;
 
-	scatterFieldAccess(scatterFieldAccess&&) = default;
+	scatteredFieldAccess(scatteredFieldAccess&&) = default;
 
-	scatterFieldAccess& operator=(const scatterFieldAccess&) = default;
+	scatteredFieldAccess& operator=(const scatteredFieldAccess&) = default;
 
-	scatterFieldAccess& operator=(scatterFieldAccess&&) = default;
+	scatteredFieldAccess& operator=(scatteredFieldAccess&&) = default;
 
-	~scatterFieldAccess() = default;
+	~scatteredFieldAccess() = default;
 
 	// - Methods
 
@@ -95,10 +95,35 @@ public:
 		return size_ == 0;
 	}
 
+	T getFirstCopy()const
+	{
+		T val;
+		uint32 n = indices_(0);
+		getNth(val, fieldVals_, n);
+		return val;
+	}
+
+	T getLastCopy()const 
+	{
+		T val;
+		if(empty())return val;
+		uint32 n = indices_(size()-1);
+		getNth(val, fieldVals_, n);
+		return val;
+	}
+
 };
 
-}
+template<typename T>
+	using deviceScatteredFieldAccess = 
+		scatteredFieldAccess<T, typename DefaultExecutionSpace::memory_space>;
+
+template<typename T>
+	using hostScatteredFieldAccess = 
+		scatteredFieldAccess<T,HostSpace>;
+
+} // pFlow
 
 
-#endif //__scatterFieldAccess_hpp__
+#endif //__scatteredFieldAccess_hpp__
 

src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp

@@ -19,8 +19,10 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "boundaryExit.hpp"
+#include "internalPoints.hpp"
 #include "dictionary.hpp"
 
+
 pFlow::boundaryExit::boundaryExit
 (
 	const dictionary& dict,
@@ -36,12 +38,120 @@ pFlow::boundaryExit::boundaryExit
 		dict.getValOrSet("checkForExitInterval", 1), 1);
 }
 
-bool pFlow::boundaryExit::beforeIteratoin
+bool pFlow::boundaryExit::beforeIteration
 (
 	uint32 iterNum, 
-	real t
+	real t,
+	real dt
 )
 {
+	// nothing have to be done
+	if(empty())
+	{
+		return true;
+	}
+
+	uint32 s = size();
+	deviceViewType1D<uint32> delFlags("delFlags",s+1); 
+	deviceViewType1D<uint32> keepFlags("keepFlags", s+1);
+	fill(delFlags, 0, s+1, 0u);
+	fill(keepFlags, 0, s+1, 0u);
+
+	using policy = Kokkos::RangePolicy<
+			pFlow::DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<pFlow::uint32>>;
+
+	auto points = thisPoints();
+	uint32 numDeleted = 0;
+	auto p = boundaryPlane().infPlane();
+
+	Kokkos::parallel_reduce
+	(
+		"boundaryExit::beforeIteration",
+		policy(0,s),
+		LAMBDA_HD(uint32 i, uint32& delToUpdate)
+		{
+			if(p.pointInNegativeSide(points(i)))
+			{
+				delFlags(i)=1;
+				delToUpdate++;
+			}
+			else
+			{
+				keepFlags(i) = 1;
+			}
+		}, 
+		numDeleted
+	);
+	
+	// no point is deleted
+	if(numDeleted == 0 )
+	{
+		return true;
+	}
+
+	exclusiveScan(delFlags, 0u, s+1, delFlags, 0u);
+	exclusiveScan(keepFlags, 0u, s+1, keepFlags, 0u);
+	
+	deviceViewType1D<uint32> deleteList("deleteList", numDeleted);
+
+	Kokkos::parallel_for
+	(
+		"boundaryExit::parllel_for",
+		policy(0, size()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(delFlags(i)!= delFlags(i+1))
+				deleteList(delFlags(i)) = i;
+		}	
+	);
+	Kokkos::fence();
+
+	deviceScatteredFieldAccess<uint32> deleteIndices(
+			numDeleted, 
+			deleteList,
+			indexList().deviceVectorAll());
+	
+	// tell internal to remove these points from its list 
+	if(!internal().deletePoints(deleteIndices))
+	{
+		fatalErrorInFunction<<
+		"error in deleting points from boundary "<< name()<<endl;
+		return false;
+	}
+
+	// delete these indices from your list 
+	if(numDeleted == s )
+	{
+		indexList().resize(0u);
+	}
+	else
+	{
+		uint newSize = s-numDeleted;
+		deviceViewType1D<uint32> newIndices("newIndices", newSize);
+		auto oldIndices = indexList().deviceVectorAll();
+
+		Kokkos::parallel_for(
+			"fillIndices",
+			policy(0,s),
+			LAMBDA_HD(uint32 i)
+			{
+				if(keepFlags(i)!= keepFlags(i+1))
+				{
+					newIndices(keepFlags(i)) = oldIndices(i);
+				}
+			}
+		);
+		Kokkos::fence();
+		copy(oldIndices, newIndices);
+		indexList().resize(newSize);
+	}
+
+	// notify your observers 
+
+	WARNING<<"notify observers in boundary exit "<<END_WARNING;
+	
 	return true;
 }
 

src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.hpp

@@ -61,7 +61,7 @@ public:
 		dictionary
 	);
 
-	bool beforeIteratoin(uint32 iterNum, real t) override;
+	bool beforeIteration(uint32 iterNum, real t, real dt) override;
 
 	bool iterate(uint32 iterNum, real t) override;
 

src/phasicFlow/structuredData/boundaries/boundaryList.cpp

@@ -19,8 +19,8 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "boundaryList.hpp"
-#include "internalPoints.hpp"
+#include "pointStructure.hpp"
-#include "simulationDomain.hpp"
+
 
 bool pFlow::boundaryList::resetLists()
 {
@@ -41,10 +41,10 @@ bool pFlow::boundaryList::updateLists()
 	dist[5] = boundary(5).neighborLength();
 
 	
-	internal_.updateFlag(
+	pStruct_.updateFlag(
-			simDomain_.thisDomain(),
+			pStruct_.simDomain().thisDomain(),
 			dist);
-	const auto& maskD = internal_.activePointsMaskDevice();
+	const auto& maskD = pStruct_.activePointsMaskDevice();
 	boundary(0).setSize( maskD.leftSize() );
 	boundary(1).setSize( maskD.rightSize() );
 	boundary(2).setSize( maskD.bottomSize() );
@@ -52,7 +52,7 @@ bool pFlow::boundaryList::updateLists()
 	boundary(4).setSize( maskD.rearSize() );
 	boundary(5).setSize( maskD.frontSize() );
 
-	internal_.fillNeighborsLists(
+	pStruct_.fillNeighborsLists(
 		boundary(0).indexList().deviceVectorAll(),
 		boundary(1).indexList().deviceVectorAll(),
 		boundary(2).indexList().deviceVectorAll(),
@@ -65,17 +65,16 @@ bool pFlow::boundaryList::updateLists()
 
 pFlow::boundaryList::boundaryList
 (
-    const simulationDomain &simD,
+    pointStructure& pStruct
-    internalPoints &internal
 )
 : 
-	ListPtr<boundaryBase>(simD.sizeOfBoundaries()),
+	ListPtr<boundaryBase>(pStruct.simDomain().sizeOfBoundaries()),
-	internal_(internal),
+	pStruct_(pStruct),
-	simDomain_(simD),
 	timeControl_
 	(
-		simDomain_.subDict("boundaries"), 
+		pStruct.simDomain().subDict("boundaries"), 
-		"update"
+		"update",
+		pStruct.currentTime()
 	)
 {}
 
@@ -92,16 +91,16 @@ bool pFlow::boundaryList::setLists()
 {
 	if(listSet_)return true;
 
-	for(auto i=0; i<simDomain_.sizeOfBoundaries();i++)
+	for(auto i=0; i<pStruct_.simDomain().sizeOfBoundaries();i++)
 	{
 		this->set
 		(
 			i,
 			boundaryBase::create
 			(
-				simDomain_.boundaryDict(i),
+				pStruct_.simDomain().boundaryDict(i),
-				simDomain_.boundaryPlane(i),
+				pStruct_.simDomain().boundaryPlane(i),
-				internal_
+				pStruct_
 			)
 		);
 	}
@@ -109,3 +108,54 @@ bool pFlow::boundaryList::setLists()
 
 	return true;
 }
+
+bool pFlow::boundaryList::beforeIteration
+(
+	uint32 iter, 
+	real t, 
+	real dt
+)
+{
+	// it is time to update lists 
+	if(timeControl_.timeEvent(iter, t, dt))
+	{
+		if(!updateLists())
+		{
+			fatalErrorInFunction;
+			return false;
+		}
+
+		WARNING<<"Maybe notification about the update list "<<END_WARNING;
+	}
+
+	for(auto& bdry:*this)
+	{
+		if( !bdry->beforeIteration(iter, t, dt))
+		{
+			fatalErrorInFunction<<
+			"Error in beforeIteration in boundary "<<bdry->name()<<endl;
+			return false;
+		}
+	}
+    return true;
+}
+
+bool pFlow::boundaryList::iterate
+(
+	uint32 iter, 
+	real t, 
+	real dt
+)
+{
+    return true;
+}
+
+bool pFlow::boundaryList::afterIteration
+(
+	uint32 iter, 
+	real t, 
+	real dt
+)
+{
+    return true;
+}

src/phasicFlow/structuredData/boundaries/boundaryList.hpp

@@ -30,8 +30,7 @@ Licence:
 namespace pFlow
 {
 
-class simulationDomain;
+class pointStructure;
-class internalPoints;
 
 class boundaryList
 :
@@ -41,13 +40,11 @@ class boundaryList
 protected:
 
 	//// - data members
-        internalPoints& 		internal_;
+		pointStructure& 		pStruct_;
-		
-        const simulationDomain& simDomain_;
 
 		baseTimeControl 		timeControl_;
 
-        bool                    listSet_ = false;
+		bool                    listSet_ = false;
 
 		bool resetLists();
 
@@ -62,9 +59,7 @@ public:
 	
 
 	//// - Constructors
-	boundaryList(
+	boundaryList(pointStructure& pStruct);
-		const simulationDomain& simD, 
-		internalPoints& internal);
 
 
   	~boundaryList() = default;
@@ -73,8 +68,6 @@ public:
 	/// the time intervals
 	bool updateLists(uint32 iter, real t, real dt);
 
-	
-
 	bool setLists();
 
 	auto& boundary(size_t i)
@@ -84,8 +77,20 @@ public:
 
 	const auto& boundary(size_t i)const
 	{
-	return ListPtr<boundaryBase>::operator[](i);
+		return ListPtr<boundaryBase>::operator[](i);
 	}
+
+	const baseTimeControl& timeControl()const
+	{
+		return timeControl_;
+	}
+
+	bool beforeIteration(uint32 iter, real t, real dt);
+
+	bool iterate(uint32 iter, real t, real dt);
+
+	bool afterIteration(uint32 iter, real t, real dt);
+	
 };
 
 } // pFlow

src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.cpp

@@ -30,10 +30,11 @@ pFlow::boundaryNone::boundaryNone
 	boundaryBase(dict, bplane, internal)
 {}
 
-bool pFlow::boundaryNone::beforeIteratoin
+bool pFlow::boundaryNone::beforeIteration
 (
 	uint32 iterNum, 
-	real t
+	real t,
+	real dt
 )
 {
 	return true;

src/phasicFlow/structuredData/boundaries/boundaryNone/boundaryNone.hpp

@@ -50,7 +50,7 @@ public:
 		dictionary
 	);
 
-	bool beforeIteratoin(uint32 iterNum, real t) override;
+	bool beforeIteration(uint32 iterNum, real t, real dt) override;
 
 	bool iterate(uint32 iterNum, real t) override;
 

src/phasicFlow/structuredData/plane/plane.hpp

@@ -122,6 +122,11 @@ public:
 	// return the parallel plane to this plane 
 	plane parallelPlane(real distance)const;
 
+	infinitePlane infPlane()const
+	{
+		return infinitePlane(normal(), d());
+	}
+
 
 	static 
 	bool validPlane4(

src/phasicFlow/structuredData/pointStructure/internalPoints.cpp

@@ -124,7 +124,21 @@ typename pFlow::internalPoints::PointsTypeHost
     return aPoints;
 }
 
-
+bool pFlow::internalPoints::deletePoints
+(
+	scatteredFieldAccess<uint32, memory_space> delPoints
+)
+{
+	if(!pFlagsD_.deletePoints(delPoints))
+	{
+		fatalErrorInFunction<<
+		"Error in deleting points from internal points"<<endl;
+		return false;
+	}
+	WARNING<<"Notify the observersin in internalPoints"<<END_WARNING;
+	pFlagSync_ = false;
+    return true;
+}
 
 FUNCTION_H
 pFlow::uint32 pFlow::internalPoints::updateFlag
@@ -133,6 +147,7 @@ pFlow::uint32 pFlow::internalPoints::updateFlag
 	const std::array<real,6>& dist
 )
 {
+	pFlagSync_ = false;
 	return pFlagsD_.markPointRegions
 	(
 		dm,
@@ -143,9 +158,7 @@ pFlow::uint32 pFlow::internalPoints::updateFlag
 		dist[3],
 		dist[4],
 		dist[5]
-	);
+	);	
-
-	pFlagSync_ = false;
 }
 
 void pFlow::internalPoints::fillNeighborsLists
@@ -166,7 +179,6 @@ void pFlow::internalPoints::fillNeighborsLists
 		rearList,
 		frontList);
 		
-	pFlagSync_ = false;
 }
 
 FUNCTION_H

src/phasicFlow/structuredData/pointStructure/internalPoints.hpp

@@ -166,6 +166,11 @@ public:
 			return pFlagsD_.activeRange();
 		}
 
+		///@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(

src/phasicFlow/structuredData/pointStructure/pointFlag.hpp

@@ -23,6 +23,7 @@ Licence:
 
 #include "phasicFlowKokkos.hpp"
 #include "domain.hpp"
+#include "scatteredFieldAccess.hpp"
 
 namespace pFlow
 {
@@ -260,6 +261,9 @@ public:
 		ViewType1D<realx3, memory_space> points);
 
 	
+	bool deletePoints(
+		scatteredFieldAccess<uint32, memory_space> points);
+	
 	/// @brief mark points based on their position in the domain.
 	/// This should be the first method to be called when updating
 	/// boundaries (step 1 of 2).

src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp

@@ -148,6 +148,55 @@ pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::scanPointFlag()
 	return numActive;
 }*/
 
+template<typename ExecutionSpace>
+bool pFlow::pointFlag<ExecutionSpace>::deletePoints
+(
+	scatteredFieldAccess<uint32, memory_space>  points
+)
+{
+	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,  
+		Kokkos::IndexType<uint32>>;
+
+	uint32 numDeleted = 0;
+	Kokkos::parallel_reduce
+	(
+		"pointFlagKernels::deletePoints",
+		policy(0u, points.size()),
+		CLASS_LAMBDA_HD(uint32 i, uint32& valDelUpdate)
+		{
+			uint32 n = points(i);
+			if(isActive(n))
+			{
+				valDelUpdate++;
+				flags_[n] = Flag::DELETED;
+			}
+		},
+		numDeleted
+	);	
+
+	if(numDeleted >= numActive_)
+	{
+		activeRange_ = {0, 0};
+		numDeleted == numActive_;
+	}
+		
+	numActive_ 	= numActive_ - numDeleted;
+	isAllActive_ = 
+		(activeRange_.numElements() == numActive_)&& numActive_>0;
+	
+	return true;
+}
+
 
 template<typename ExecutionSpace>
 pFlow::uint32 pFlow::pointFlag<ExecutionSpace>::markPointRegions

src/phasicFlow/structuredData/pointStructure/pointStructure.cpp

@@ -114,17 +114,17 @@ pFlow::pointStructure::pointStructure
 	),
 	boundaries_
 	(
-		simulationDomain_(),
+        *this
-		*this
 	)
 {
-    REPORT(0)<< "Reading point structure from "<<
+    REPORT(0)<< "Reading "<< Green_Text("point structure")<<
-    IOobject::path()<<END_REPORT;
+    " from "<<IOobject::path()<<END_REPORT;
 
     if( !IOobject::readObject() )
     {
         fatalErrorInFunction<<
         "Error in reading from file "<<IOobject::path()<<endl;
+        fatalExit;
     }
 }
 
@@ -153,7 +153,6 @@ pFlow::pointStructure::pointStructure(
 	),
 	boundaries_
 	(
-		simulationDomain_(),
 		*this
 	)
 {