add an additional P2P: P2POuter which compute the P2P between neighbors...

add an additional P2P: P2POuter which compute the P2P between neighbors (possibly in mutual way) but not the inner P2P in the target leaf

Src/Adaptive/FAdaptChebSymKernel.hpp

@@ -55,7 +55,7 @@ class FTreeCoordinate;
  */
 
 template<class FReal, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
-class FAdaptiveChebSymKernel : FChebSymKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
+class FAdaptiveChebSymKernel : public FChebSymKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
         , public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
     //
     typedef FChebSymKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>	KernelBaseClass;

Src/Adaptive/FAdaptUnifKernel.hpp

@@ -56,7 +56,7 @@ class FTreeCoordinate;
  */
 
 template<class FReal, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
-class FAdaptiveUnifKernel : FUnifKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
+class FAdaptiveUnifKernel : public FUnifKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
         , public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
     //
     typedef FUnifKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>	KernelBaseClass;

Src/Adaptive/FAdaptiveKernelWrapper.hpp

@@ -347,6 +347,13 @@ public:
         kernel.P2P(inLeafPosition, targets, sources, directNeighborsParticles, positions, size);
 	}
 
+    /** This is a normal P2P */
+    void P2POuter(const FTreeCoordinate& inLeafPosition,
+            ContainerClass* const FRestrict targets,
+            ContainerClass* const directNeighborsParticles[], const int positions[], const int size)  override {
+        kernel.P2POuter(inLeafPosition, targets, directNeighborsParticles, positions, size);
+    }
+
 	/** This is a normal P2P */
 	void P2PRemote(const FTreeCoordinate& inLeafPosition,
 			ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,

Src/BalanceTree/FChebSymCostKernel.hpp

@@ -322,6 +322,34 @@ public:
         cell->addNearCost(tmpCost);
         countP2P++;
     }
+
+    void P2POuter(const FTreeCoordinate& LeafCellCoordinate, // needed for periodic boundary conditions
+             ContainerClass* const FRestrict TargetParticles,
+             ContainerClass* const NeighborSourceParticles[],
+             const int positions[],
+             const int size) override {
+        FSize tmpCost = 0;
+        FSize tgtPartCount = TargetParticles->getNbParticles();
+
+        {
+            for (int idx=0; idx < size && positions[idx]<=13; ++idx)
+            {
+                    tmpCost +=
+                        countFlopsP2Pmutual()
+                        * tgtPartCount
+                        * NeighborSourceParticles[idx]->getNbParticles();
+                }
+        }
+
+        flopsP2P += tmpCost;
+
+        CellClass* cell = _tree->getCell(
+            LeafCellCoordinate.getMortonIndex(_treeHeight - 1),
+            _treeHeight - 1);
+
+        cell->addNearCost(tmpCost);
+        countP2P++;
+    }
 };
 
 

Src/BalanceTree/FDumbCostKernel.hpp

@@ -193,6 +193,22 @@ public:
         FSize tmpCost = srcPartCount * tgtPartCount;
 
 
+        CellClass* cell = _tree->getCell(
+            LeafCellCoordinate.getMortonIndex(_treeHeight - 1),
+            _treeHeight - 1);
+
+        flopsP2P += tmpCost;
+        cell->addNearCost(tmpCost);
+        countP2P++;
+    }
+
+    void P2POuter(const FTreeCoordinate& LeafCellCoordinate, // needed for periodic boundary conditions
+             ContainerClass* const FRestrict TargetParticles,
+             ContainerClass* const /*NeighborSourceParticles*/[],
+             const int /*positions*/[],
+             const int size) override {
+        FSize tmpCost = 0;
+
         CellClass* cell = _tree->getCell(
             LeafCellCoordinate.getMortonIndex(_treeHeight - 1),
             _treeHeight - 1);

Src/Components/FAbstractKernels.hpp

@@ -27,9 +27,10 @@
 * @brief This class defines what any kernel has to implement.
 *
 * Please notice that P2PRemote is optional and should be implemented in case of
-* MPI usage.
+* MPI usage and certain algorithm.
 * It is better to inherit from this class even if it is not obligatory thanks to
-* the templates. But inheriting will force your class to have the correct parameters.
+* the templates. But inheriting will force your class to have the correct parameters
+* especially of you use the override keyword.
 *
 * You can find an example of implementation in FBasicKernels.
 */
@@ -41,44 +42,55 @@ public:
     }
 
     /**
-        * P2M
-        * particles to multipole
-        * @param pole the multipole to fill using the particles
-        * @param particles the particles from the same spacial boxe
-        */
+    * P2M
+    * particles to multipole
+    * @param pole the multipole to fill using the particles
+    * @param particles the particles from the same spacial boxe
+    */
     virtual void P2M(CellClass* const pole, const ContainerClass* const particles) = 0;
 
     /**
-        * M2M
-        * Multipole to multipole
-        * @param pole the father (the boxe that contains other ones)
-        * @param child the boxe to take values from
-        * @param the current computation level
-        * the child array has a size of 8 elements (address if exists or 0 otherwise).
-        * You must test if a pointer is 0 to know if an element exists inside this array
-        */
+    * M2M
+    * Multipole to multipole
+    * @param pole the father (the boxe that contains other ones)
+    * @param child the boxe to take values from
+    * @param the current computation level
+    * the child array has a size of 8 elements (address if exists or 0 otherwise).
+    * You must test if a pointer is 0 to know if an element exists inside this array.
+    * The order of the M2M are driven by the Morton indexing (the children appear in the
+    * array in increasing Morton index).
+    * In binary: xyz, so 100 means 1 in x, 0 in y, 0 in z.
+    */
     virtual void M2M(CellClass* const FRestrict pole, const CellClass*const FRestrict *const FRestrict child, const int inLevel) = 0;
 
     /**
-        * M2L
-        * Multipole to local
-        * @param local the element to fill using distant neighbors
-        * @param distantNeighbors is an array containing fathers's direct neighbors's child - direct neigbors (max 189)
-        * @param neighborPositions the relative position of the neighbor (between O and 342)
-        * @param size the number of neighbors
-        * @param inLevel the current level of the computation
-        *
-        * @code // the code was:
-        * @code  for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
-        * @code      if(distantNeighbors[idxNeigh]){
-        * @code           ...
-        * @code      }
-        * @code  }
-        * @code  // Now it must be :
-        * @code  for(int idxExistingNeigh = 0 ; idxExistingNeigh < size ; ++idxExistingNeigh){
-        * @code      const int idxNeigh = neighborPositions[idxExistingNeigh]
-        * @code      distantNeighbors[idxExistingNeigh]...
-        * @code  }
+    * M2L
+    * Multipole to local
+    * @param local the element to fill using distant neighbors
+    * @param distantNeighbors is an array containing fathers's direct neighbors's child - direct neigbors
+    *        (max 189 in normal M2L and 189+26 with extended separation criteria)
+    * @param neighborPositions the relative position of the neighbor (between O and 342)
+    * @param size the number of neighbors
+    * @param inLevel the current level of the computation
+    *
+    * The relative position is given by the formula:
+    * (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3
+    * with xdiff, ydiff and zdiff from -3 to 2 or from -2 to 3.
+    *
+    * The way of accessing the interaction list has been changed in a previous ScalFMM:
+    * @code // the code was:
+    * @code  for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
+    * @code      // Test if a cell exists
+    * @code      if(distantNeighbors[idxNeigh]){
+    * @code           ...
+    * @code      }
+    * @code  }
+    * @code  // Now it must be :
+    * @code  for(int idxExistingNeigh = 0 ; idxExistingNeigh < size ; ++idxExistingNeigh){
+    * @code      const int idxNeigh = neighborPositions[idxExistingNeigh]
+    * @code      distantNeighbors[idxExistingNeigh]...
+    * @code  }
+    * It may have negligable extra cost in dense FMM but is clearly benefits for Sparse FMM.
         */
     virtual void M2L(CellClass* const FRestrict local, const CellClass* distantNeighbors[],
                      const int neighborPositions[],
@@ -86,19 +98,21 @@ public:
 
 
     /** This method is used to bypass needFinishedM2LEvent method a each level
-       *  during the transferPass.  If you have to use the  finishedLevelM2L then you
-       *  have to  inherit it and return true rather than false.
-       *
-       * @return false
-       */
-     constexpr static bool NeedFinishedM2LEvent(){
-    	 	 return false;
-     }
-     /** This method can be optionally inherited
+    *  during the transferPass.  If you have to use the  finishedLevelM2L then you
+    *  have to  inherit it and return true rather than false.
+    * But it will imply extra dependencies and even not be supported by some algorithms.
+    *
+    * @return false
+    */
+    constexpr static bool NeedFinishedM2LEvent(){
+        return false;
+    }
+
+    /** This method can be optionally inherited
       * It is called at the end of each computation level during the M2L pass
       * @param level the ending level
       */
-    void finishedLevelM2L(const int /*level*/){
+    virtual void finishedLevelM2L(const int /*level*/){
     }
 
     /**
@@ -108,8 +122,10 @@ public:
         * @param child the child to downward values (child may have already been impacted by M2L)
         * @param inLevel the current level of computation
         * the child array has a size of 8 elements (address if exists or 0 otherwise).
-        * Children are ordering in the morton index way.
-        * You must test if a pointer is 0 to know if an element exists inside this array
+        * You must test if a pointer is 0 to know if an element exists inside this array.
+        * The order of the M2M are driven by the Morton indexing (the children appear in the
+        * array in increasing Morton index).
+        * In binary: xyz, so 100 means 1 in x, 0 in y, 0 in z.
         */
     virtual void L2L(const CellClass* const FRestrict local, CellClass* FRestrict * const FRestrict child, const int inLevel) = 0;
 
@@ -124,36 +140,84 @@ public:
     /**
         * P2P
         * Particles to particles
-        * @param inLeafPosition tree coordinate of the leaf
+        * This functions should compute the inner interactions inside the target leaf,
+        * and the interaction between the target and the list of neighbors.
+        * This neighbor computation should be the same as P2POuter and thus,
+        * internally this function may call P2POuter directly.
+        *
+        * The neighborPositions contains a value from 0 to 26 included which represent the relative
+        * position of the neighbor given by:
+        * (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1
+        * with idxX, idxY, and idxZ from -1 to 1.
+        *
+        *
+        * @param inLeafPosition tree coordinate of the leaf (number of boxes in x, y and z)
         * @param targets current boxe targets particles
         * @param sources current boxe sources particles (can be == to targets)
-        * @param directNeighborsParticles the particles from direct neighbors (this is an array of list)  max length = 26
-        * @param neighborPositions the relative position of the neighbors (between O and 25) max length = 26
+        * @param directNeighborsParticles the particles from direct neighbors (max length = 26)
+        * @param neighborPositions the relative position of the neighbors (between O and 25, max length = 26)
         * @param size the number of direct neighbors
         */
     virtual void P2P(const FTreeCoordinate& inLeafPosition,
-             ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-             ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
-             const int size) = 0;
+                     ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
+                     ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
+                     const int size) = 0;
 
     /**
         * P2P
         * Particles to particles
-        * @param inLeafPosition tree coordinate of the leaf
+        * This functions should compute the interaction between the target and the list of neighbors.
+        * This neighbor computation should be the same as P2P and thus,
+        * this function may be called by P2P.
+        *
+        * The neighborPositions contains a value from 0 to 26 included which represent the relative
+        * position of the neighbor given by:
+        * (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1
+        * with idxX, idxY, and idxZ from -1 to 1.
+        *
+        *
+        * @param inLeafPosition tree coordinate of the leaf (number of boxes in x, y and z)
         * @param targets current boxe targets particles
-        * @param sources current boxe sources particles (can be == to targets)
-        * @param directNeighborsParticles the particles from direct neighbors (this is an array of list) max length = 26
-        * @param neighborPositions the relative position of the neighbors (between 0 and 25)  max length = 26
+        * @param directNeighborsParticles the particles from direct neighbors (max length = 26)
+        * @param neighborPositions the relative position of the neighbors (between O and 25, max length = 26)
         * @param size the number of direct neighbors
-        *
-        * This method is called by the MPI algorithm with leaves from other hosts.
-        * Notice that directNeighborsParticles will be destroyed once all P2P remote have been
-        * performed.
         */
-  virtual void P2PRemote(const FTreeCoordinate& /*inLeafPosition*/,
-			 ContainerClass* const FRestrict /*targets*/, const ContainerClass* const FRestrict /*sources*/,
-             ContainerClass* const /*directNeighborsParticles*/[],
-             const int /*neighborPositions*/[], const int /*size*/) {
+    virtual void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+                          ContainerClass* const FRestrict /*targets*/,
+                          ContainerClass* const /*directNeighborsParticles*/[], const int /*neighborPositions*/[],
+                          const int /*size*/)  = 0;
+
+    /**
+    * P2P
+    * Particles to particles
+    * This functions should compute the interaction between the target and the list of neighbors.
+    *
+    * There are no interest to compute mutual interaction involve in this function.
+    * It is mainly called by the MPI algorithm with leaves from other hosts and thus
+    * modifying the neighbors do not modify the original data:
+    * directNeighborsParticles will be destroyed once all P2P remote have been performed.
+    *
+    * The neighborPositions contains a value from 0 to 26 included which represent the relative
+    * position of the neighbor given by:
+    * (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1
+    * with idxX, idxY, and idxZ from -1 to 1.
+    *
+    * The sources is given but the inner P2P must not be computed here.
+    * It is given just in case the P2PRemote need it, but the interaction between
+    * targets and sources should be done in the P2P function.
+    *
+    * @param inLeafPosition tree coordinate of the leaf (number of boxes in x, y and z)
+    * @param targets current boxe targets particles
+    * @param sources current boxe sources particles (can be == to targets)
+    * @param directNeighborsParticles the particles from direct neighbors (max length = 26)
+    * @param neighborPositions the relative position of the neighbors (between O and 25, max length = 26)
+    * @param size the number of direct neighbors
+    *
+    */
+    virtual void P2PRemote(const FTreeCoordinate& /*inLeafPosition*/,
+                           ContainerClass* const FRestrict /*targets*/, const ContainerClass* const FRestrict /*sources*/,
+                           ContainerClass* const /*directNeighborsParticles*/[],
+                           const int /*neighborPositions*/[], const int /*size*/) {
         FLOG( FLog::Controller.write("Warning, P2P remote is used but not implemented!").write(FLog::Flush) );
     }
 

Src/Components/FBasicKernels.hpp

@@ -68,6 +68,14 @@ public:
 
     }
 
+    /** Do nothing */
+    virtual void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+                          ContainerClass* const FRestrict /*targets*/,
+                          ContainerClass* const /*directNeighborsParticles*/[], const int /*neighborPositions*/[],
+                          const int /*size*/) override {
+
+    }
+
     /** Do nothing */
     virtual void P2PRemote(const FTreeCoordinate& /*treeCoord*/,
                            ContainerClass* const FRestrict /*targetParticles*/, const ContainerClass* const FRestrict /*sourceParticles*/,

Src/Components/FTestKernels.hpp

@@ -111,6 +111,22 @@ public:
         }
     }
 
+    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+             ContainerClass* const FRestrict targets,
+             ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
+             const int inSize) override {
+        long long int inc = 0;
+
+        for(int idx = 0 ; idx < inSize ; ++idx){
+                inc += directNeighborsParticles[idx]->getNbParticles();
+        }
+
+        long long int*const particlesAttributes = targets->getDataDown();
+        for(FSize idxPart = 0 ; idxPart < targets->getNbParticles() ; ++idxPart){
+            particlesAttributes[idxPart] += inc;
+        }
+    }
+
     /** After Downward */
     void P2PRemote(const FTreeCoordinate& ,
                    ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,

Src/GroupTree/Core/FGroupSeqAlgorithm.hpp

@@ -385,12 +385,12 @@ protected:
                                 FAssertLF((*iterParticles)->getLeafMortonIndex(outsideInteractions[outInterIdx].insideIdxInBlock) == outsideInteractions[outInterIdx].insideIndex);
 
                                 ParticleContainerClass* ptrLeaf = &interParticles;
-                                kernels->P2PRemote( FTreeCoordinate(outsideInteractions[outInterIdx].insideIndex, tree->getHeight()-1),
-                                                    &particles, &particles , &ptrLeaf, &outsideInteractions[outInterIdx].outPosition, 1);
+                                kernels->P2POuter( FTreeCoordinate(outsideInteractions[outInterIdx].insideIndex, tree->getHeight()-1),
+                                                    &particles , &ptrLeaf, &outsideInteractions[outInterIdx].outPosition, 1);
                                 const int otherPosition = getOppositeNeighIndex(outsideInteractions[outInterIdx].outPosition);
                                 ptrLeaf = &particles;
-                                kernels->P2PRemote( FTreeCoordinate(outsideInteractions[outInterIdx].outIndex, tree->getHeight()-1),
-                                                    &interParticles, &interParticles , &ptrLeaf, &otherPosition, 1);
+                                kernels->P2POuter( FTreeCoordinate(outsideInteractions[outInterIdx].outIndex, tree->getHeight()-1),
+                                                    &interParticles , &ptrLeaf, &otherPosition, 1);
                             }
                         }
                     }

Src/GroupTree/Core/FGroupTaskAlgorithm.hpp

@@ -653,12 +653,12 @@ protected:
                                 FAssertLF(containers->getLeafMortonIndex((*outsideInteractions)[outInterIdx].insideIdxInBlock) == (*outsideInteractions)[outInterIdx].insideIndex);                              
 
                                 ParticleContainerClass* ptrLeaf = &interParticles;
-                                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, tree->getHeight()-1),
-                                                    &particles, &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
+                                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, tree->getHeight()-1),
+                                                    &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
                                 const int otherPosition = getOppositeNeighIndex((*outsideInteractions)[outInterIdx].outPosition);
                                 ptrLeaf = &particles;
-                                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, tree->getHeight()-1),
-                                                    &interParticles, &interParticles , &ptrLeaf, &otherPosition, 1);
+                                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, tree->getHeight()-1),
+                                                    &interParticles , &ptrLeaf, &otherPosition, 1);
                             }
                         }
                     }

Src/GroupTree/Core/FGroupTaskDepAlgorithm.hpp

@@ -811,12 +811,12 @@ protected:
                                 FAssertLF(containers->getLeafMortonIndex((*outsideInteractions)[outInterIdx].insideIdxInBlock) == (*outsideInteractions)[outInterIdx].insideIndex);
 
                                 ParticleContainerClass* ptrLeaf = &interParticles;
-                                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, tree->getHeight()-1),
-                                                    &particles, &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
+                                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, tree->getHeight()-1),
+                                                    &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
                                 const int otherPosition = getOppositeNeighIndex((*outsideInteractions)[outInterIdx].outPosition);
                                 ptrLeaf = &particles;
-                                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, tree->getHeight()-1),
-                                                    &interParticles, &interParticles , &ptrLeaf, &otherPosition, 1);
+                                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, tree->getHeight()-1),
+                                                    &interParticles , &ptrLeaf, &otherPosition, 1);
                             }
                         }
                     }

Src/GroupTree/StarPUUtils/FStarPUCpuWrapper.hpp

@@ -489,12 +489,12 @@ public:
                 FAssertLF(containers->getLeafMortonIndex((*outsideInteractions)[outInterIdx].insideIdxInBlock) == (*outsideInteractions)[outInterIdx].insideIndex);
 
                 ParticleContainerClass* ptrLeaf = &interParticles;
-                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, treeHeight-1),
-                                    &particles, &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
+                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].insideIndex, treeHeight-1),
+                                    &particles , &ptrLeaf, &(*outsideInteractions)[outInterIdx].outPosition, 1);
                 const int otherPosition = getOppositeNeighIndex((*outsideInteractions)[outInterIdx].outPosition);
                 ptrLeaf = &particles;
-                kernel->P2PRemote( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, treeHeight-1),
-                                    &interParticles, &interParticles , &ptrLeaf, &otherPosition, 1);
+                kernel->P2POuter( FTreeCoordinate((*outsideInteractions)[outInterIdx].outIndex, treeHeight-1),
+                                    &interParticles , &ptrLeaf, &otherPosition, 1);
             }
         }
     }

Src/Kernels/Chebyshev/FAbstractChebKernel.hpp

@@ -156,6 +156,11 @@ public:
             const int SourcePositions[],
 		   const int /* size */) = 0;
 
+  virtual void P2POuter(const FTreeCoordinate& inLeafPosition,
+           ContainerClass* const FRestrict targets,
+           ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
+           const int size) = 0;
+
 
   virtual void P2PRemote(const FTreeCoordinate& /*inPosition*/,
 			 ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,

Src/Kernels/Chebyshev/FChebFlopsSymKernel.hpp

@@ -205,10 +205,10 @@ public:
         unsigned int flops = 0;
 		// count how ofter each of the 16 interactions is used
 		memset(countExp, 0, sizeof(int) * 343);
-        for (unsigned int idx=0; idx<size; ++idx)
+        for (int idx=0; idx<size; ++idx)
             countExp[SymHandler->pindices[positions[idx]]]++;
 		// multiply (mat-mat-mul)
-        for (unsigned int pidx=0; pidx<343; ++pidx)
+        for (int pidx=0; pidx<343; ++pidx)
 			if (countExp[pidx])
 				flops += countFlopsM2L(countExp[pidx], SymHandler->LowRank[pidx]) + countExp[pidx]*nnodes;
 		flopsM2L += flops;
@@ -254,15 +254,24 @@ public:
 	{
 		if (TargetParticles != SourceParticles) {
             flopsP2P += countFlopsP2P() * TargetParticles->getNbParticles() * SourceParticles->getNbParticles();
-            for (unsigned int idx=0; idx<size; ++idx)
+            for (int idx=0; idx<size; ++idx)
                     flopsP2P += countFlopsP2P() * TargetParticles->getNbParticles() * NeighborSourceParticles[idx]->getNbParticles();
 		} else {
             flopsP2P += countFlopsP2Pmutual() * ((TargetParticles->getNbParticles()*TargetParticles->getNbParticles()+TargetParticles->getNbParticles()) / 2);
-            for (unsigned int idx=0; idx < size && neighborPositions[idx]<=13; ++idx)
+            for (int idx=0; idx < size && neighborPositions[idx]<=13; ++idx)
                     flopsP2P += countFlopsP2Pmutual() * TargetParticles->getNbParticles() * NeighborSourceParticles[idx]->getNbParticles();
 		}
 	}
 
+    void P2POuter(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
+                     ContainerClass* const FRestrict TargetParticles,
+                     ContainerClass* const NeighborSourceParticles[],
+                     const int neighborPositions[],
+                     const int size) override
+    {
+        for (int idx=0; idx < size && neighborPositions[idx]<=13; ++idx)
+                flopsP2P += countFlopsP2Pmutual() * TargetParticles->getNbParticles() * NeighborSourceParticles[idx]->getNbParticles();
+    }
 
 };
 

Src/Kernels/Chebyshev/FChebKernel.hpp

@@ -206,10 +206,18 @@ public:
 
     }
 
-    void P2P(const FTreeCoordinate& /*inPosition*/,
+    void P2P(const FTreeCoordinate& inPosition,
              ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
              ContainerClass* const inNeighbors[], const int neighborPositions[],
              const int inSize) override {
+        P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
+        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PInner(inTargets,MatrixKernel);
+    }
+
+    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+             ContainerClass* const FRestrict inTargets,
+             ContainerClass* const inNeighbors[], const int neighborPositions[],
+             const int inSize) override {
         int nbNeighborsToCompute = 0;
         while(nbNeighborsToCompute < inSize
               && neighborPositions[nbNeighborsToCompute] < 14){
@@ -218,7 +226,6 @@ public:
         DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2P(inTargets,inNeighbors,nbNeighborsToCompute,MatrixKernel);
     }
 
-
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[], const int /*neighborPositions*/[],

Src/Kernels/Chebyshev/FChebSymKernel.hpp

@@ -447,10 +447,18 @@ public:
 
     }
 
-    void P2P(const FTreeCoordinate& /*inPosition*/,
+    void P2P(const FTreeCoordinate& inPosition,
              ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
              ContainerClass* const inNeighbors[], const int neighborPositions[],
              const int inSize) override {
+        P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
+        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PInner(inTargets,MatrixKernel);
+    }
+
+    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+             ContainerClass* const FRestrict inTargets,
+             ContainerClass* const inNeighbors[], const int neighborPositions[],
+             const int inSize) override {
         int nbNeighborsToCompute = 0;
         while(nbNeighborsToCompute < inSize
               && neighborPositions[nbNeighborsToCompute] < 14){

Src/Kernels/Chebyshev/FChebTensorialKernel.hpp

@@ -219,10 +219,18 @@ public:
         }
     }
 
-    void P2P(const FTreeCoordinate& /*inPosition*/,
+    void P2P(const FTreeCoordinate& inPosition,
              ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
              ContainerClass* const inNeighbors[], const int neighborPositions[],
              const int inSize) override {
+        P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
+        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PInner(inTargets,MatrixKernel);
+    }
+
+    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
+             ContainerClass* const FRestrict inTargets,
+             ContainerClass* const inNeighbors[], const int neighborPositions[],
+             const int inSize) override {
         int nbNeighborsToCompute = 0;
         while(nbNeighborsToCompute < inSize
               && neighborPositions[nbNeighborsToCompute] < 14){

Src/Kernels/Interpolation/FInterpP2PKernels.hpp

@@ -20,6 +20,12 @@ struct DirectInteractionComputer
       FP2P::FullMutualKIJ<FReal, ContainerClass, MatrixKernelClass>(TargetParticles,NeighborSourceParticles,inSize,MatrixKernel);
   }
 
+  template <typename ContainerClass, typename MatrixKernelClass>
+  static void P2PInner( ContainerClass* const FRestrict TargetParticles,
+                   const MatrixKernelClass *const MatrixKernel){