From 3e7bcd0a860b670f3367f1fa9aa4f66ab81e9c13 Mon Sep 17 00:00:00 2001
From: berenger-bramas <berenger-bramas@2616d619-271b-44dc-8df4-d4a8f33a7222>
Date: Mon, 13 Feb 2012 10:51:37 +0000
Subject: [PATCH] Change P2P prototype.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@377 2616d619-271b-44dc-8df4-d4a8f33a7222
---
 Src/Components/FAbstractKernels.hpp          |  13 +--
 Src/Components/FBasicKernels.hpp             |   9 +-
 Src/Components/FTestKernels.hpp              |  29 ++---
 Src/Containers/FOctree.hpp                   |   8 +-
 Src/Core/FFmmAlgorithm.hpp                   |   4 +-
 Src/Core/FFmmAlgorithmPeriodic.hpp           |   4 +-
 Src/Core/FFmmAlgorithmTask.hpp               |  18 +--
 Src/Core/FFmmAlgorithmThread.hpp             |   4 +-
 Src/Core/FFmmAlgorithmThreadProc.hpp         |  30 +++--
 Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp |  32 +++---
 Src/Core/FFmmAlgorithmThreadTsm.hpp          |   5 +-
 Src/Core/FFmmAlgorithmTsm.hpp                |   5 +-
 Src/Kernels/FAbstractSphericalKernel.hpp     | 112 ++++---------------
 Tests/testFmmDemonstration.cpp               |  13 ++-
 Tests/testStatsTree.cpp                      |   2 +-
 15 files changed, 120 insertions(+), 168 deletions(-)

diff --git a/Src/Components/FAbstractKernels.hpp b/Src/Components/FAbstractKernels.hpp
index c9fcc19cd..6b0f7095d 100644
--- a/Src/Components/FAbstractKernels.hpp
+++ b/Src/Components/FAbstractKernels.hpp
@@ -82,26 +82,25 @@ public:
         * P2P
         * Particles to particles
         * @param targets current boxe targets particles
-        * @param sources current boxe sources particles
         * @param directNeighborsParticles the particles from direct neighbors (this is an array of list)
         * @param size the number of direct neighbors (the size of the array directNeighborsParticles)
         */
-    virtual void P2P(ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-                     const ContainerClass* const directNeighborsParticles[26], const int size) = 0;
+    virtual void P2P(const FTreeCoordinate& inLeafPosition,
+             ContainerClass* const FRestrict targets,
+             ContainerClass* const directNeighborsParticles[27], const int size) = 0;
 
     /**
         * P2P
         * Particles to particles
-        * @param inCurrentLeafIndex
+        * @param inLeafPosition
         * @param targets current boxe targets particles
         * @param sources current boxe sources particles
         * @param directNeighborsParticles the particles from direct neighbors (this is an array of list)
-        * @param inNeighborsIndex the indexes of neighbors
         * @param size the number of direct neighbors (the size of the array directNeighborsParticles)
         */
-    virtual void P2P(const MortonIndex inCurrentLeafIndex,
+    virtual void P2P(const FTreeCoordinate& inLeafPosition,
              ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-             ContainerClass* const directNeighborsParticles[26], const int size) = 0;
+             ContainerClass* const directNeighborsParticles[27], const int size) = 0;
 
 };
 
diff --git a/Src/Components/FBasicKernels.hpp b/Src/Components/FBasicKernels.hpp
index c9aca4e99..47b3ae8f6 100644
--- a/Src/Components/FBasicKernels.hpp
+++ b/Src/Components/FBasicKernels.hpp
@@ -59,15 +59,16 @@ public:
     }
 
     /** Do nothing */
-    virtual void P2P(ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
-                     const ContainerClass* const [26], const int ) {
+    virtual void P2P(const FTreeCoordinate& ,
+                     ContainerClass* const FRestrict ,
+                     ContainerClass* const [27], const int ) {
 
     }
 
     /** Do nothing */
-    virtual void P2P(const MortonIndex ,
+    virtual void P2P(const FTreeCoordinate& ,
                      ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
-                     ContainerClass* const [26], const int ){
+                     ContainerClass* const [27], const int ){
 
     }
 
diff --git a/Src/Components/FTestKernels.hpp b/Src/Components/FTestKernels.hpp
index c8f8d500a..02273c14c 100644
--- a/Src/Components/FTestKernels.hpp
+++ b/Src/Components/FTestKernels.hpp
@@ -90,16 +90,17 @@ public:
     }
 
     /** After Downward */
-    void P2P(ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-             const ContainerClass* const directNeighborsParticles[26], const int size) {
+    void P2P(const FTreeCoordinate& ,
+                 ContainerClass* const FRestrict targets,
+                 ContainerClass* const directNeighborsParticles[27], const int ){
 
         // Each particles targeted is impacted by the particles sources
-        long long int inc = sources->getSize();
-        if(targets == sources){
-            inc -= 1;
-        }
-        for(int idx = 0 ; idx < size ; ++idx){
-            inc += directNeighborsParticles[idx]->getSize();
+        long long int inc = targets->getSize() - 1;
+
+        for(int idx = 0 ; idx < 27 ; ++idx){
+            if( directNeighborsParticles[idx] ){
+                inc += directNeighborsParticles[idx]->getSize();
+            }
         }
 
         typename ContainerClass::BasicIterator iter(*targets);
@@ -112,17 +113,19 @@ public:
 
 
     /** After Downward */
-    void P2P(const MortonIndex ,
-             ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-             ContainerClass* const directNeighborsParticles[26], const int size) {
+    void P2P(const FTreeCoordinate& ,
+                 ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
+                 ContainerClass* const directNeighborsParticles[27], const int ){
 
         // Each particles targeted is impacted by the particles sources
         long long int inc = sources->getSize();
         if(targets == sources){
             inc -= 1;
         }
-        for(int idx = 0 ; idx < size ; ++idx){
-            inc += directNeighborsParticles[idx]->getSize();
+        for(int idx = 0 ; idx < 27 ; ++idx){
+            if( directNeighborsParticles[idx] ){
+                inc += directNeighborsParticles[idx]->getSize();
+            }
         }
 
         typename ContainerClass::BasicIterator iter(*targets);
diff --git a/Src/Containers/FOctree.hpp b/Src/Containers/FOctree.hpp
index 2ffbed481..074441f25 100644
--- a/Src/Containers/FOctree.hpp
+++ b/Src/Containers/FOctree.hpp
@@ -933,8 +933,8 @@ public:
           * @param inLevel the level of the element
           * @return the number of neighbors
           */
-    int getLeafsNeighbors(ContainerClass* inNeighbors[26], const FTreeCoordinate& center, const int inLevel){
-        memset( inNeighbors, 0 , 26 * sizeof(ContainerClass*));
+    int getLeafsNeighbors(ContainerClass* inNeighbors[27], const FTreeCoordinate& center, const int inLevel){
+        memset( inNeighbors, 0 , 27 * sizeof(ContainerClass*));
         const int limite = FMath::pow2(inLevel);
 
         int idxNeighbors = 0;
@@ -974,8 +974,8 @@ public:
           * @param inLevel the level of the element
           * @return the number of neighbors
           */
-    int getPeriodicLeafsNeighbors(ContainerClass* inNeighbors[26], const FTreeCoordinate& center, const int inLevel){
-        memset(inNeighbors , 0 , 26 * sizeof(ContainerClass*));
+    int getPeriodicLeafsNeighbors(ContainerClass* inNeighbors[27], const FTreeCoordinate& center, const int inLevel){
+        memset(inNeighbors , 0 , 27 * sizeof(ContainerClass*));
         const int limite = FMath::pow2(inLevel);
 
         int idxNeighbors = 0;
diff --git a/Src/Core/FFmmAlgorithm.hpp b/Src/Core/FFmmAlgorithm.hpp
index 6176814ec..4216b4dbb 100644
--- a/Src/Core/FFmmAlgorithm.hpp
+++ b/Src/Core/FFmmAlgorithm.hpp
@@ -233,7 +233,7 @@ private:
         typename OctreeClass::Iterator octreeIterator(tree);
         octreeIterator.gotoBottomLeft();
         // There is a maximum of 26 neighbors
-        ContainerClass* neighbors[26];
+        ContainerClass* neighbors[27];
         // for each leafs
         do{
             FDEBUG(computationCounterL2P.tic());
@@ -242,7 +242,7 @@ private:
             // need the current particles and neighbors particles
             const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),heightMinusOne);
             FDEBUG(computationCounterP2P.tic());
-            kernels->P2P(octreeIterator.getCurrentGlobalIndex(),octreeIterator.getCurrentListTargets(), octreeIterator.getCurrentListSrc() , neighbors, counter);
+            kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(), neighbors, counter);
             FDEBUG(computationCounterP2P.tac());
         } while(octreeIterator.moveRight());
 
diff --git a/Src/Core/FFmmAlgorithmPeriodic.hpp b/Src/Core/FFmmAlgorithmPeriodic.hpp
index 1a005e55a..980fab517 100644
--- a/Src/Core/FFmmAlgorithmPeriodic.hpp
+++ b/Src/Core/FFmmAlgorithmPeriodic.hpp
@@ -231,7 +231,7 @@ private:
         typename OctreeClass::Iterator octreeIterator(tree);
         octreeIterator.gotoBottomLeft();
         // There is a maximum of 26 neighbors
-        ContainerClass* neighbors[26];
+        ContainerClass* neighbors[27];
         // for each leafs
         do{
             FDEBUG(computationCounterL2P.tic());
@@ -240,7 +240,7 @@ private:
             // need the current particles and neighbors particles
             const int counter = tree->getPeriodicLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),heightMinusOne);
             FDEBUG(computationCounterP2P.tic());
-            kernels->P2P(octreeIterator.getCurrentGlobalIndex(),octreeIterator.getCurrentListTargets(), octreeIterator.getCurrentListSrc() , neighbors, counter);
+            kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(), neighbors, counter);
             FDEBUG(computationCounterP2P.tac());
         } while(octreeIterator.moveRight());
 
diff --git a/Src/Core/FFmmAlgorithmTask.hpp b/Src/Core/FFmmAlgorithmTask.hpp
index 7a6bac76d..61d2d6e7d 100644
--- a/Src/Core/FFmmAlgorithmTask.hpp
+++ b/Src/Core/FFmmAlgorithmTask.hpp
@@ -277,13 +277,13 @@ private:
 
         const int heightMinusOne = OctreeHeight - 1;
 
-#pragma omp parallel
+        #pragma omp parallel
         {
             KernelClass*const myThreadkernels = kernels[omp_get_thread_num()];
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
 
-#pragma omp single nowait
+            #pragma omp single nowait
             {
                 const int SizeShape = 3*3*3;
                 FVector<typename OctreeClass::Iterator> shapes[SizeShape];
@@ -293,7 +293,7 @@ private:
 
                 // for each leafs
                 do{
-#pragma omp task
+                    #pragma omp task
                     {
                         myThreadkernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
                     }
@@ -302,7 +302,7 @@ private:
                     const int shapePosition = (coord.getX()%3)*9 + (coord.getY()%3)*3 + (coord.getZ()%3);
 
                     if( shapePosition == 0){
-#pragma omp task
+                        #pragma omp task
                         {
                             // need the current particles and neighbors particles
                             const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),heightMinusOne);
@@ -315,20 +315,20 @@ private:
 
                 } while(octreeIterator.moveRight());
 
-#pragma omp taskwait
+                #pragma omp taskwait
 
                 for( int idxShape = 1 ; idxShape < SizeShape ; ++idxShape){
                     int iterLeaf = shapes[idxShape].getSize();
                     while( iterLeaf-- ){
                         typename OctreeClass::Iterator toWork = shapes[idxShape][iterLeaf];
-#pragma omp task
+                        #pragma omp task
                         {
                             const int counter = tree->getLeafsNeighbors(neighbors, toWork.getCurrentGlobalCoordinate(),heightMinusOne);
-                            myThreadkernels->P2P(toWork.getCurrentGlobalIndex(),toWork.getCurrentListTargets(), toWork.getCurrentListSrc() , neighbors, counter);
+                            myThreadkernels->P2P(toWork.getCurrentGlobalCoordinate(), toWork.getCurrentListTargets(), neighbors, counter);
                         }
                     }
 
-#pragma omp taskwait
+                    #pragma omp taskwait
                 }
             }
         }
diff --git a/Src/Core/FFmmAlgorithmThread.hpp b/Src/Core/FFmmAlgorithmThread.hpp
index be47737e5..2c34c7e14 100644
--- a/Src/Core/FFmmAlgorithmThread.hpp
+++ b/Src/Core/FFmmAlgorithmThread.hpp
@@ -371,7 +371,7 @@ private:
 
             KernelClass& myThreadkernels = (*kernels[omp_get_thread_num()]);
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
             int previous = 0;
 
             for(int idxShape = 0 ; idxShape < SizeShape ; ++idxShape){
@@ -384,7 +384,7 @@ private:
                     // need the current particles and neighbors particles
                     FDEBUG(if(!omp_get_thread_num()) computationCounterP2P.tic());
                     const int counter = tree->getLeafsNeighbors(neighbors, currentIter.cell->getCoordinate(), LeafIndex);
-                    myThreadkernels.P2P(currentIter.index, currentIter.targets, currentIter.sources , neighbors, counter);
+                    myThreadkernels.P2P(currentIter.cell->getCoordinate(), currentIter.targets, neighbors, counter);
                     FDEBUG(if(!omp_get_thread_num()) computationCounterP2P.tac());
                 }
 
diff --git a/Src/Core/FFmmAlgorithmThreadProc.hpp b/Src/Core/FFmmAlgorithmThreadProc.hpp
index 43abacf61..b7a9b1a23 100644
--- a/Src/Core/FFmmAlgorithmThreadProc.hpp
+++ b/Src/Core/FFmmAlgorithmThreadProc.hpp
@@ -953,6 +953,7 @@ private:
             int alreadySent[nbProcess];
 
             MortonIndex indexesNeighbors[26];
+            int uselessIndexArray[26];
 
             for(int idxLeaf = 0 ; idxLeaf < this->numberOfLeafs ; ++idxLeaf){
                 FTreeCoordinate center;
@@ -961,7 +962,7 @@ private:
                 memset(alreadySent, 0, sizeof(int) * nbProcess);
                 bool needOther = false;
 
-                const int neighCount = getNeighborsIndexes(iterArray[idxLeaf].getCurrentGlobalIndex(), limite, indexesNeighbors);
+                const int neighCount = getNeighborsIndexes(iterArray[idxLeaf].getCurrentGlobalIndex(), limite, indexesNeighbors, uselessIndexArray);
 
                 for(int idxNeigh = 0 ; idxNeigh < neighCount ; ++idxNeigh){
                     if(indexesNeighbors[idxNeigh] < intervals[idProcess].min || intervals[idProcess].max < indexesNeighbors[idxNeigh]){
@@ -1118,7 +1119,7 @@ private:
         {
             KernelClass& myThreadkernels = (*kernels[omp_get_thread_num()]);
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
             int previous = 0;
 
             for(int idxShape = 0 ; idxShape < SizeShape ; ++idxShape){
@@ -1132,7 +1133,7 @@ private:
 
                         // need the current particles and neighbors particles
                         const int counter = tree->getLeafsNeighbors(neighbors, currentIter.cell->getCoordinate(), LeafIndex);
-                        myThreadkernels.P2P( currentIter.index,currentIter.targets, currentIter.sources , neighbors, counter);
+                        myThreadkernels.P2P( currentIter.cell->getCoordinate(),currentIter.targets, neighbors, counter);
                     }
                 }
 
@@ -1181,11 +1182,11 @@ private:
         FTRACE( FTrace::FRegion regionOtherTrace("Compute P2P Other", __FUNCTION__ , __FILE__ , __LINE__) );
         FDEBUG( computation2Counter.tic() );
 
-#pragma omp parallel
+        #pragma omp parallel
         {
             KernelClass& myThreadkernels = (*kernels[omp_get_thread_num()]);
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
             int previous = 0;
             // Box limite
             const int limite = 1 << (this->OctreeHeight - 1);
@@ -1204,20 +1205,23 @@ private:
                         int counter = tree->getLeafsNeighbors(neighbors, currentIter.cell->getCoordinate(), LeafIndex);
 
                         // Take possible data
-                        MortonIndex indexesNeighbors[26];
-                        const int neighCount = getNeighborsIndexes(currentIter.index, limite, indexesNeighbors);
+                        MortonIndex indexesNeighbors[27];
+                        int indexArray[26];
+                        const int nbNeigh = getNeighborsIndexes(currentIter.index, limite, indexesNeighbors, indexArray);
+
+                        memset( neighbors, 0, sizeof(ContainerClass*) * 27);
 
-                        for(int idxNeigh = 0 ; idxNeigh < neighCount ; ++idxNeigh){
+                        for(int idxNeigh = 0 ; idxNeigh < nbNeigh ; ++idxNeigh){
                             if(indexesNeighbors[idxNeigh] < intervals[idProcess].min || intervals[idProcess].max < indexesNeighbors[idxNeigh]){
                                 ContainerClass*const hypotheticNeighbor = otherP2Ptree.getLeafSrc(indexesNeighbors[idxNeigh]);
                                 if(hypotheticNeighbor){
-                                    neighbors[counter] = hypotheticNeighbor;
+                                    neighbors[ indexArray[idxNeigh] ] = hypotheticNeighbor;
                                     ++counter;
                                 }
                             }
                         }
 
-                        myThreadkernels.P2P( currentIter.index,currentIter.targets, currentIter.sources , neighbors, counter);
+                        myThreadkernels.P2P( currentIter.cell->getCoordinate(), currentIter.targets, neighbors, counter);
                     }
                 }
 
@@ -1245,7 +1249,7 @@ private:
     }
 
 
-    int getNeighborsIndexes(const MortonIndex centerIndex, const int limite, MortonIndex indexes[26]) const{
+    int getNeighborsIndexes(const MortonIndex centerIndex, const int limite, MortonIndex indexes[26], int indexInArray[26]) const{
         FTreeCoordinate center;
         center.setPositionFromMorton(centerIndex, OctreeHeight - 1);
         int idxNeig = 0;
@@ -1262,7 +1266,9 @@ private:
                     // if we are not on the current cell
                     if( idxX || idxY || idxZ ){
                         const FTreeCoordinate other(center.getX() + idxX,center.getY() + idxY,center.getZ() + idxZ);
-                        indexes[idxNeig++] = other.getMortonIndex(this->OctreeHeight - 1);
+                        indexes[ idxNeig ] = other.getMortonIndex(this->OctreeHeight - 1);
+                        indexInArray[ idxNeig ] = ((idxX+1)*3 + (idxY+1)) * 3 + (idxZ+1);
+                        ++idxNeig;
                     }
                 }
             }
diff --git a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
index 6a78f0159..616e66314 100644
--- a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
+++ b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
@@ -1031,6 +1031,7 @@ private:
             int alreadySent[nbProcess];
 
             MortonIndex indexesNeighbors[26];
+            int uselessIndexInArray[26];
 
             for(int idxLeaf = 0 ; idxLeaf < this->numberOfLeafs ; ++idxLeaf){
                 FTreeCoordinate center;
@@ -1039,9 +1040,9 @@ private:
                 memset(alreadySent, 0, sizeof(int) * nbProcess);
                 bool needOther = false;
 
-                getNeighborsIndexes(iterArray[idxLeaf].getCurrentGlobalCoordinate(), limite, indexesNeighbors);
+                const int nbNeigh = getNeighborsIndexes(iterArray[idxLeaf].getCurrentGlobalCoordinate(), limite, indexesNeighbors, uselessIndexInArray);
 
-                for(int idxNeigh = 0 ; idxNeigh < 26 ; ++idxNeigh){
+                for(int idxNeigh = 0 ; idxNeigh < nbNeigh ; ++idxNeigh){
                     if(indexesNeighbors[idxNeigh] < intervals[idProcess].min || intervals[idProcess].max < indexesNeighbors[idxNeigh]){
                         needOther = true;
 
@@ -1195,7 +1196,7 @@ private:
         {
             KernelClass& myThreadkernels = (*kernels[omp_get_thread_num()]);
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
             int previous = 0;
 
             for(int idxShape = 0 ; idxShape < SizeShape ; ++idxShape){
@@ -1209,7 +1210,7 @@ private:
 
                         // need the current particles and neighbors particles
                         const int counter = tree->getPeriodicLeafsNeighbors(neighbors, currentIter.cell->getCoordinate(), LeafIndex);
-                        myThreadkernels.P2P( currentIter.index,currentIter.targets, currentIter.sources , neighbors, counter);
+                        myThreadkernels.P2P( currentIter.cell->getCoordinate(), currentIter.targets, neighbors, counter);
                     }
                 }
 
@@ -1262,8 +1263,9 @@ private:
         {
             KernelClass& myThreadkernels = (*kernels[omp_get_thread_num()]);
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
             MortonIndex indexesNeighbors[26];
+            int indexInArray[26];
 
             int previous = 0;
             // Box limite
@@ -1282,28 +1284,28 @@ private:
                         // need the current particles and neighbors particles
                         // int counter = tree->getLeafsNeighbors(neighbors, neighborsIndex, currentIter.index, LeafIndex);
                         int counter = 0;
-                        memset(neighbors, 0, 26 * sizeof(ContainerClass*));
+                        memset(neighbors, 0, 27 * sizeof(ContainerClass*));
                         // Take possible data
-                        getNeighborsIndexes(currentIter.cell->getCoordinate(), limite, indexesNeighbors);
+                        const int nbNeigh = getNeighborsIndexes(currentIter.cell->getCoordinate(), limite, indexesNeighbors, indexInArray);
 
-                        for(int idxNeigh = 0 ; idxNeigh < 26 ; ++idxNeigh){
+                        for(int idxNeigh = 0 ; idxNeigh < nbNeigh ; ++idxNeigh){
                             if(indexesNeighbors[idxNeigh] < intervals[idProcess].min || intervals[idProcess].max < indexesNeighbors[idxNeigh]){
                                 ContainerClass*const hypotheticNeighbor = otherP2Ptree.getLeafSrc(indexesNeighbors[idxNeigh]);
                                 if(hypotheticNeighbor){
-                                    neighbors[idxNeigh] = hypotheticNeighbor;
+                                    neighbors[ indexInArray[idxNeigh] ] = hypotheticNeighbor;
                                     ++counter;
                                 }
                             }
                             else{
                                 ContainerClass*const hypotheticNeighbor = tree->getLeafSrc(indexesNeighbors[idxNeigh]);
                                 if(hypotheticNeighbor){
-                                    neighbors[idxNeigh] = hypotheticNeighbor;
+                                    neighbors[ indexInArray[idxNeigh] ] = hypotheticNeighbor;
                                     ++counter;
                                 }
                             }
                         }
 
-                        myThreadkernels.P2P( currentIter.index,currentIter.targets, currentIter.sources , neighbors, counter);
+                        myThreadkernels.P2P( currentIter.cell->getCoordinate(), currentIter.targets, neighbors, counter);
                     }
                 }
 
@@ -1332,7 +1334,7 @@ private:
 
 
 
-    void getNeighborsIndexes(const FTreeCoordinate& center, const int limite, MortonIndex indexes[26]) const{
+    int getNeighborsIndexes(const FTreeCoordinate& center, const int limite, MortonIndex indexes[26], int indexInArray[26]) const{
         int idxNeig = 0;
         // We test all cells around
         for(int idxX = -1 ; idxX <= 1 ; ++idxX){
@@ -1352,11 +1354,13 @@ private:
                         else if( limite <= other.getZ() ) other.setZ( other.getZ() - limite );
 
                         indexes[idxNeig] = other.getMortonIndex(this->OctreeHeight - 1);
-                    }
-                    ++idxNeig;
+                        indexInArray[idxNeig] = ((idxX+1)*3  + (idxY+1))*3 + (idxZ+1);
+                        ++idxNeig;
+                    }                    
                 }
             }
         }
+        return idxNeig;
     }
 
 
diff --git a/Src/Core/FFmmAlgorithmThreadTsm.hpp b/Src/Core/FFmmAlgorithmThreadTsm.hpp
index 44d8c75fa..baca291e9 100644
--- a/Src/Core/FFmmAlgorithmThreadTsm.hpp
+++ b/Src/Core/FFmmAlgorithmThreadTsm.hpp
@@ -343,14 +343,15 @@ public:
         {
             KernelClass * const myThreadkernels = kernels[omp_get_thread_num()];
             // There is a maximum of 26 neighbors
-            ContainerClass* neighbors[26];
+            ContainerClass* neighbors[27];
 
             #pragma omp for schedule(dynamic)
             for(int idxLeafs = 0 ; idxLeafs < numberOfLeafs ; ++idxLeafs){
                 myThreadkernels->L2P(iterArray[idxLeafs].getCurrentCell(), iterArray[idxLeafs].getCurrentListTargets());
                 // need the current particles and neighbors particles
                 const int counter = tree->getLeafsNeighbors(neighbors, iterArray[idxLeafs].getCurrentGlobalCoordinate(),heightMinusOne);
-                myThreadkernels->P2P( iterArray[idxLeafs].getCurrentListTargets(), iterArray[idxLeafs].getCurrentListSrc() , neighbors, counter);
+                myThreadkernels->P2P( iterArray[idxLeafs].getCurrentGlobalCoordinate(), iterArray[idxLeafs].getCurrentListTargets(),
+                                      iterArray[idxLeafs].getCurrentListSrc() , neighbors, counter);
             }
         }
         FDEBUG(computationCounter.tac());
diff --git a/Src/Core/FFmmAlgorithmTsm.hpp b/Src/Core/FFmmAlgorithmTsm.hpp
index 6198e120b..aefef1e72 100644
--- a/Src/Core/FFmmAlgorithmTsm.hpp
+++ b/Src/Core/FFmmAlgorithmTsm.hpp
@@ -274,14 +274,15 @@ public:
         typename OctreeClass::Iterator octreeIterator(tree);
         octreeIterator.gotoBottomLeft();
         // There is a maximum of 26 neighbors
-        ContainerClass* neighbors[26];
+        ContainerClass* neighbors[27];
         // for each leafs
         do{
             FDEBUG(computationCounter.tic());
             kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
             // need the current particles and neighbors particles
             const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), heightMinusOne);
-            kernels->P2P( octreeIterator.getCurrentListTargets(), octreeIterator.getCurrentListSrc() , neighbors, counter);
+            kernels->P2P( octreeIterator.getCurrentGlobalCoordinate(), octreeIterator.getCurrentListTargets(),
+                          octreeIterator.getCurrentListSrc() , neighbors, counter);
             FDEBUG(computationCounter.tac());
             FDEBUG(totalComputation += computationCounter.elapsed());
         } while(octreeIterator.moveRight());
diff --git a/Src/Kernels/FAbstractSphericalKernel.hpp b/Src/Kernels/FAbstractSphericalKernel.hpp
index 6827ff78c..7b3b1e078 100644
--- a/Src/Kernels/FAbstractSphericalKernel.hpp
+++ b/Src/Kernels/FAbstractSphericalKernel.hpp
@@ -197,8 +197,9 @@ public:
       * then it computes the sources/targets inteactions between this leaf and the
       * neighbors.
       */
-    void P2P(ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
-             const ContainerClass* const directNeighbors[26], const int size) {
+    void P2P(const FTreeCoordinate& inLeafPosition,
+                  ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
+                  ContainerClass* const directNeighborsParticles[27], const int ){
 
         { // Compute interaction in this leaf
             typename ContainerClass::BasicIterator iterTarget(*targets);
@@ -209,7 +210,6 @@ public:
                 // For all the source particles in the same leaf
                 typename ContainerClass::ConstBasicIterator iterSameBox(*sources);
                 while( iterSameBox.hasNotFinished() ){
-                    //(&iterSameBox.data() != &iterTarget.data())
                     directInteraction(&target, iterSameBox.data());
                     iterSameBox.gotoNext();
                 }
@@ -220,66 +220,16 @@ public:
         }
         { // Compute interactions with other leaves
             // For all the neigbors leaves
-            for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
-                // For all particles in current leaf
-                typename ContainerClass::BasicIterator iterTarget(*targets);
-                while( iterTarget.hasNotFinished() ){
-                    ParticleClass target( iterTarget.data() );
-                    // For all the particles in the other leaf
-                    typename ContainerClass::ConstBasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
-                    while( iterSource.hasNotFinished() ){
-                        directInteraction(&target, iterSource.data());
-                        iterSource.gotoNext();
-                    }
-                    // Set data and progress
-                    iterTarget.setData(target);
-                    iterTarget.gotoNext();
-                }
-            }
-        }
-    }
-
-    /** This P2P has to be used when target == sources
-      * It will proceed a direct interation >> mutual
-      *
-      * It takes all the particles from the current leaf,
-      * then it computes the interactions in this leaf,
-      * then it computes the  inteactions between this leaf and the
-      * neighbors.
-      */
-    void P2P(const MortonIndex inCurrentIndex,
-             ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
-             ContainerClass* const directNeighbors[26], const int size){
-        { // Compute interaction in this leaf
-            typename ContainerClass::BasicIterator iterTarget(*targets);
-            while( iterTarget.hasNotFinished() ){
-                // We copy the target particle to work with a particle in the heap
-                ParticleClass target( iterTarget.data() );
-
-                // For all particles after the current one
-                typename ContainerClass::BasicIterator iterSameBox = iterTarget;
-                iterSameBox.gotoNext();
-                while( iterSameBox.hasNotFinished() ){
-                    directInteractionMutual(&target, &iterSameBox.data());
-                    iterSameBox.gotoNext();
-                }
-                // Set data and progress
-                iterTarget.setData(target);
-                iterTarget.gotoNext();
-            }
-        }
-        { // Compute interactions with other leaves
-            // For all the neigbors leaves
-            for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
-                if(inCurrentIndex < 0 /*inNeighborsIndex[idxDirectNeighbors] TODO */ ){
+            for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < 27 ; ++idxDirectNeighbors){
+                if( directNeighborsParticles[idxDirectNeighbors] ){
                     // For all particles in current leaf
                     typename ContainerClass::BasicIterator iterTarget(*targets);
                     while( iterTarget.hasNotFinished() ){
                         ParticleClass target( iterTarget.data() );
                         // For all the particles in the other leaf
-                        typename ContainerClass::BasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
+                        typename ContainerClass::ConstBasicIterator iterSource(*directNeighborsParticles[idxDirectNeighbors]);
                         while( iterSource.hasNotFinished() ){
-                            directInteractionMutual(&target, &iterSource.data());
+                            directInteraction(&target, iterSource.data());
                             iterSource.gotoNext();
                         }
                         // Set data and progress
@@ -291,16 +241,17 @@ public:
         }
     }
 
-    ///////////////////////////////////////////////////////////////////////////////
-    //                                  Periodic
-    ///////////////////////////////////////////////////////////////////////////////
-
-
-    /** After Downward */
-    /*void P2P(const MortonIndex inCurrentIndex,
-             ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,//sources
-             ContainerClass* const directNeighbors[26], const FTreeCoordinate neighborsRelativeOffset[26], const int size) {
-
+    /** This P2P has to be used when target == sources
+      * It will proceed a direct interation >> mutual
+      *
+      * It takes all the particles from the current leaf,
+      * then it computes the interactions in this leaf,
+      * then it computes the  inteactions between this leaf and the
+      * neighbors.
+      */
+    void P2P(const FTreeCoordinate& inLeafPosition,
+             ContainerClass* const FRestrict targets,
+             ContainerClass* const directNeighborsParticles[27], const int ){
         { // Compute interaction in this leaf
             typename ContainerClass::BasicIterator iterTarget(*targets);
             while( iterTarget.hasNotFinished() ){
@@ -321,32 +272,14 @@ public:
         }
         { // Compute interactions with other leaves
             // For all the neigbors leaves
-            for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
-                // This box is not a real neighbor
-                if(neighborsRelativeOffset[idxDirectNeighbors].getX() || neighborsRelativeOffset[idxDirectNeighbors].getY()
-                        || neighborsRelativeOffset[idxDirectNeighbors].getZ() ){
-                    typename ContainerClass::BasicIterator iterTarget(*targets);
-                    while( iterTarget.hasNotFinished() ){
-                        ParticleClass target( iterTarget.data() );
-                        // For all the particles in the other leaf
-                        typename ContainerClass::BasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
-                        while( iterSource.hasNotFinished() ){
-                            directInteractionOffset(&target, iterSource.data(), neighborsRelativeOffset[idxDirectNeighbors]);
-                            iterSource.gotoNext();
-                        }
-                        // Set data and progress
-                        iterTarget.setData(target);
-                        iterTarget.gotoNext();
-                    }
-                }
-                // This is a real neighbor, we do as usual
-                else if(inCurrentIndex < neighborsRelativeOffset[idxDirectNeighbors].getMortonIndex(treeHeight) ){
+            for(int idxDirectNeighbors = 0 ; idxDirectNeighbors <= 13 ; ++idxDirectNeighbors){
+                if( directNeighborsParticles[idxDirectNeighbors] ){
                     // For all particles in current leaf
                     typename ContainerClass::BasicIterator iterTarget(*targets);
                     while( iterTarget.hasNotFinished() ){
                         ParticleClass target( iterTarget.data() );
                         // For all the particles in the other leaf
-                        typename ContainerClass::BasicIterator iterSource(*directNeighbors[idxDirectNeighbors]);
+                        typename ContainerClass::BasicIterator iterSource(*directNeighborsParticles[idxDirectNeighbors]);
                         while( iterSource.hasNotFinished() ){
                             directInteractionMutual(&target, &iterSource.data());
                             iterSource.gotoNext();
@@ -358,7 +291,8 @@ public:
                 }
             }
         }
-    }*/
+    }
+
 
     ///////////////////////////////////////////////////////////////////////////////
     //                                  Computation
diff --git a/Tests/testFmmDemonstration.cpp b/Tests/testFmmDemonstration.cpp
index 8e760f5a5..fff3729d8 100644
--- a/Tests/testFmmDemonstration.cpp
+++ b/Tests/testFmmDemonstration.cpp
@@ -75,7 +75,7 @@ public:
 // My kernel actually does nothing except showing how to retreive data from an
 // array from the indexes vector giving by the leaf in the P2M
 template< class ParticleClass, class CellClass, class ContainerClass>
-class MyKernel /*: public FAbstractKernels<ParticleClass,CellClass,ContainerClass> */{
+class MyKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass>{
     FBasicParticle*const realParticles;
 public:
     MyKernel(FBasicParticle*const inRealParticles): realParticles(inRealParticles) {
@@ -101,13 +101,16 @@ public:
     void L2P(const CellClass* const , ContainerClass* const ){
     }
 
-    void P2P(ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
-                     const ContainerClass* const [26], const int ) {
+    void P2P(const FTreeCoordinate& ,
+                     ContainerClass* const FRestrict ,
+                     ContainerClass* const [27], const int ) {
+
     }
 
-    void P2P(const MortonIndex ,
+    void P2P(const FTreeCoordinate& ,
                      ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
-                     ContainerClass* const [26], const int ){
+                     ContainerClass* const [27], const int ){
+
     }
 };
 
diff --git a/Tests/testStatsTree.cpp b/Tests/testStatsTree.cpp
index 208331792..d7a0d9aec 100644
--- a/Tests/testStatsTree.cpp
+++ b/Tests/testStatsTree.cpp
@@ -151,7 +151,7 @@ int main(int argc, char ** argv){
             octreeIterator.gotoBottomLeft();
 
             do{
-                ContainerClass* neighbors[26];
+                ContainerClass* neighbors[27];
                 // need the current particles and neighbors particles
                 averageNeighbors += FReal(tree.getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),NbLevels-1));
                 ++nbLeafs;
-- 
GitLab