contact searcch for Hrchl parallel

src/Interaction/contactSearch/methods/NBS.H

@@ -57,6 +57,35 @@ public:
 
 	struct TagFindPairs{};
 
+	class cellIterator
+	{
+	private:
+		ViewType3D<int32, memory_space> 		head_;
+
+		ViewType1D<int32, memory_space> 		next_;
+
+	public:
+
+		cellIterator(ViewType3D<int32, memory_space> head, ViewType1D<int32, memory_space> next)
+		:
+			head_(head),
+			next_(next)
+		{}
+
+		INLINE_FUNCTION_HD
+		Cells cellsSize()const {
+			return Cells(head_.extent(0), head_.extent(1), head_.extent(2));}
+
+		INLINE_FUNCTION_HD
+		int32 start(indexType i, indexType j, indexType k)const {
+			return head_(i,j,k); }
+
+		INLINE_FUNCTION_HD
+		int32 getNext(int32 n)const {
+			if(n<0) return n;
+			return next_(n); }
+	};
+
 protected:
 
 	int32 		capacity_		= 1;
@@ -256,6 +285,11 @@ public:
 		return next_;
 	}
 
+	cellIterator getCellIterator()const
+	{
+		return cellIterator(head_, next_);
+	}
+
 	// - Perform the broad search to find pairs
 	//	 with force = true, perform broad search regardless of 
 	//   updateFrequency_ value
@@ -480,84 +514,6 @@ public:
 		return true;
 	}
 
-
-	/*template <typename PairsContainer, typename teamMemberType>
-	INLINE_FUNCTION_HD
-	int32 addPointsInBoxToList( 
-		const teamMemberType& teamMember,
-		IdType id,
-		const iBox<IndexType>& triBox,
-		const PairsContainer& pairs)const
-	{
-		int32 getFull = 0;
-
-		auto bExtent = boxExtent(triBox);
-		int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
-
-		const auto head = head_;
-		const auto next = next_;
-
-		// perform a loop over all cells in the triBox
-		Kokkos::parallel_reduce(
-			Kokkos::TeamThreadRange(teamMember,numCellBox),
-			[=](int32 linIndex, int32& valToUpdate){
-
-				CellType cell;
-				indexToCell(linIndex, triBox, cell);
-
-				int32 n = head(cell.x(),cell.y(),cell.z());
-				
-				while( n>-1)
-				{
-					
-					// id is wall id the pair is (particle id, wall id)
-					if( pairs.insert(static_cast<IdType>(n), id) < 0 )
-						valToUpdate++;
-					n = next(n);
-				}
-
-			},
-			getFull
-			);
-
-		return getFull;
-
-	}*/
-
-	template <typename PairsContainer>
-	INLINE_FUNCTION_HD
-	int32 addPointsInBoxToListModified( 
-		IdType id,
-		const iBox<IndexType>& triBox,
-		const PairsContainer& pairs)const
-	{
-		int32 getFull = 0;
-
-		auto bExtent = boxExtent(triBox);
-		int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
-
-		const auto head = head_;
-		const auto next = next_;
-
-		for(int32 linIndex=0; linIndex<numCellBox; linIndex++)
-		{
-				CellType cell;
-				indexToCell(linIndex, triBox, cell);
-
-				int32 n = head_(cell.x(),cell.y(),cell.z());
-				while( n>-1)
-				{
-					// id is wall id the pair is (particle id, wall id)
-					if( pairs.insert(static_cast<IdType>(n), id) < 0 )
-						getFull++;
-					n = next_(n);
-				}
-		}
-		
-		return getFull;
-
-	}
-
 };
 
 }

src/Interaction/contactSearch/wallMappings/cellsSimple.H

@@ -257,7 +257,7 @@ public:
 		while (getFull)
 		{
 			
-			getFull = findPairsElementRangeModified(pairs, particleMap);
+			getFull = findPairsElementRange(pairs, particleMap);
 			
 			if(getFull)
 			{
@@ -271,8 +271,6 @@ public:
 						oldCap << " to " 
 						<< pairs.capacity() <<" in cellsSimple."<<endInfo;
 				
-				
-				
 				Kokkos::fence();
 			}
 		}
@@ -280,7 +278,8 @@ public:
 		return true;
 	}
 
-	/*template<typename PairsContainer, typename particleMapType>
+
+	template<typename PairsContainer, typename particleMapType>
 	int32 findPairsElementRange(PairsContainer& pairs, particleMapType& particleMap) 
 	{
 		int32 getFull =0;
@@ -288,65 +287,51 @@ public:
 		const auto pwPairs = pairs;
 		const auto elementBox = elementBox_;
 
+		auto cellIter = particleMap.getCellIterator();
 		
 		Kokkos::parallel_reduce(
-			"cellsSimple::findPairsElementRange",
+			"cellsSimple::findPairsElementRangeModified2",
 			tpPWContactSearch(numElements_, Kokkos::AUTO),
 			LAMBDA_HD(
-				const typename tpPWContactSearch::member_type & teamMember,
+				const typename tpPWContactSearch::member_type & teamMember,	int32& valueToUpdate){
-				int32& valueToUpdate){
-				int32 i  = teamMember.league_rank();
-				printf("league size %d , league rank %d , team size %d and team rank %d \n", teamMember.league_size(), i, teamMember.team_size(), teamMember.team_rank());
-				IdType id = i;
-				const auto triBox = elementBox[i];
 				
-				valueToUpdate +=
+				const int32 iTri = teamMember.league_rank();
-					particleMap.addPointsInBoxToList(
+
-						teamMember,
+				const auto triBox = elementBox[iTri];
-						id,
+				
-						triBox,
+				int32 getFull2 = 0;
-						pwPairs
+
-						);
+				auto bExtent = boxExtent(triBox);
+				int32 numCellBox = bExtent.x()*bExtent.y()*bExtent.z();
+
+				Kokkos::parallel_reduce( 
+					Kokkos::TeamThreadRange( teamMember, numCellBox ),
+					[&] ( const int32 linIndex, int32 &innerUpdate )
+						{
+							
+							CellType cell;
+							indexToCell(linIndex, triBox, cell);
+
+							int32 n = cellIter.start(cell.x(),cell.y(),cell.z());
+						
+							while( n>-1)
+							{
+								// id is wall id the pair is (particle id, wall id)
+								if( pairs.insert(static_cast<IdType>(n), iTri) < 0 )
+									innerUpdate++;
+								n = cellIter.getNext(n);
+							}
+
+						},
+					getFull2
+				);
+
+				if ( teamMember.team_rank() == 0 ) valueToUpdate += getFull2;
 			},
 			getFull
 			);
 
 		return getFull;
-	}*/
-
-	template<typename PairsContainer, typename particleMapType>
-	int32 findPairsElementRangeModified(PairsContainer& pairs, particleMapType& particleMap)
-	{
-		Kokkos::RangePolicy<
-			Kokkos::IndexType<int32>,
-			Kokkos::Schedule<Kokkos::Dynamic>,
-			ExecutionSpace> rPolicy(0,numElements_);
-
-		int32 getFull =0;
-		
-		const auto pwPairs = pairs;
-		const auto elementBox = elementBox_;
-		//printf("we are in modified version \n");
-		Kokkos::parallel_reduce(
-			"cellsSimple::findPairsElementRangeModified",
-			rPolicy,
-			LAMBDA_HD(
-				int32 i,
-				int32& valueToUpdate){
-				
-				IdType id = i;
-				const auto triBox = elementBox[i];
-				
-				valueToUpdate +=
-					particleMap.addPointsInBoxToListModified(
-						id,
-						triBox,
-						pwPairs
-						);
-			},
-			getFull
-			);
-		return getFull;
 	}
 
 	INLINE_FUNCTION_HD
@@ -360,7 +345,6 @@ public:
 		realx3 minP, maxP;
 
 		this->extendBox(p1, p2, p3, cellExtent_, minP, maxP);
-		//output<< minP << " maxP "<< maxP<<endl;
 		elementBox_[i] = iBoxType(this->pointIndex(minP), this->pointIndex(maxP));
 
 	}