From b86f1ffdf1429d315a8f939bae67ff37c283c307 Mon Sep 17 00:00:00 2001
From: Berenger Bramas <Berenger.Bramas@inria.fr>
Date: Mon, 27 Apr 2015 13:39:05 +0200
Subject: [PATCH] Add a way to choose the separation criteria (2 or 1 -
 default- or 0 or -1)

---
 Src/Containers/FOctree.hpp                   | 29 ++++---------
 Src/Containers/FTreeCoordinate.hpp           |  9 ++--
 Src/Core/FFmmAlgorithm.hpp                   | 10 +++--
 Src/Core/FFmmAlgorithmPeriodic.hpp           |  8 ++--
 Src/Core/FFmmAlgorithmSectionTask.hpp        |  8 ++--
 Src/Core/FFmmAlgorithmTask.hpp               |  8 ++--
 Src/Core/FFmmAlgorithmThread.hpp             |  9 ++--
 Src/Core/FFmmAlgorithmThreadProc.hpp         | 26 ++++++++----
 Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp | 43 +++++++++++++-------
 Src/Core/FFmmAlgorithmThreadTsm.hpp          |  9 ++--
 Src/Core/FFmmAlgorithmTsm.hpp                |  9 ++--
 11 files changed, 100 insertions(+), 68 deletions(-)

diff --git a/Src/Containers/FOctree.hpp b/Src/Containers/FOctree.hpp
index 607e3b1f7..c5bd0cb7f 100644
--- a/Src/Containers/FOctree.hpp
+++ b/Src/Containers/FOctree.hpp
@@ -808,8 +808,8 @@ public:
      * @return the number of neighbors
      */
     int getInteractionNeighbors(const CellClass* inNeighbors[343],
-    const FTreeCoordinate& workingCell,
-    const int inLevel) const{
+                const FTreeCoordinate& workingCell,
+                const int inLevel, const int neighSeparation = 1) const{
         // reset
         memset(inNeighbors, 0, sizeof(CellClass*) * 343);
 
@@ -848,7 +848,7 @@ public:
                                     const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
 
                                     // Test if it is a direct neighbor
-                                    if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                                    if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                         // add to neighbors
                                         inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
                                         ++idxNeighbors;
@@ -875,8 +875,8 @@ public:
      * @return the number of neighbors
      */
     int getFullNeighborhood(const CellClass* inNeighbors[343],
-    const FTreeCoordinate& workingCell,
-    const int inLevel) const{
+                const FTreeCoordinate& workingCell,
+                const int inLevel) const{
         // reset
         memset(inNeighbors, 0, sizeof(CellClass*) * 343);
 
@@ -935,10 +935,8 @@ public:
      * @return the number of neighbors
      */
     int getPeriodicInteractionNeighbors(const CellClass* inNeighbors[343],
-    const FTreeCoordinate& workingCell,
-    const int inLevel, const int inDirection) const{
-        // TODO : REMOVE NEXT COMMENTS
-        //		std::cout << " Begin in getPeriodicInteractionNeighbors"<<std::endl;
+                const FTreeCoordinate& workingCell,
+                const int inLevel, const int inDirection, const int neighSeparation = 1) const{
 
         // Then take each child of the parent's neighbors if not in directNeighbors
         // Father coordinate
@@ -963,11 +961,7 @@ public:
             const int endY =    (TestPeriodicCondition(inDirection, DirPlusY)  || parentCell.getY() != boxLimite - 1 ?1:0);
             const int startZ =  (TestPeriodicCondition(inDirection, DirMinusZ) || parentCell.getZ() != 0 ?-1:0);
             const int endZ =    (TestPeriodicCondition(inDirection, DirPlusZ)  || parentCell.getZ() != boxLimite - 1 ?1:0);
-            // TODO : REMOVE NEXT COMMENTS
-            //			std::cout << "  -- startX " << startX << " endX "<< endX<< std::endl ;
-            //			std::cout << "  -- startY " << startY << " endX "<< endY<< std::endl ;
-            //			std::cout << "  -- startZ " << startZ << " endX "<< endZ<< std::endl ;
-            //			std::cout << "     boxLimite "<< boxLimite<<std::endl;
+
             int idxNeighbors = 0;
             // We test all cells around
             for(int idxX = startX ; idxX <= endX ; ++idxX){
@@ -1018,15 +1012,10 @@ public:
                                         const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1))         - workingCell.getZ();
 
                                         // Test if it is a direct neighbor
-                                        if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                                        if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                             // add to neighbors
-                                            // TODO : REMOVE NEXT COMMENTS
-                                            //											std::cout << "     Voisin numero "<< idxNeighbors
-                                            //													  << " indexinTab "<< (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3
-                                            //													  << "  idxXousin " << idxCousin<< std::endl;
                                             inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
                                             ++idxNeighbors;
-
                                         }
                                     }
                                 }
diff --git a/Src/Containers/FTreeCoordinate.hpp b/Src/Containers/FTreeCoordinate.hpp
index f84fc8e1e..d904cc4fe 100644
--- a/Src/Containers/FTreeCoordinate.hpp
+++ b/Src/Containers/FTreeCoordinate.hpp
@@ -342,7 +342,8 @@ public:
         return idxNeig;
     }
 
-    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[189], int inNeighborsPosition[189]) const{
+    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], int inNeighborsPosition[/*189+26+1*/216],
+                            const int neighSeparation = 1) const{
         // Then take each child of the parent's neighbors if not in directNeighbors
         // Father coordinate
         const FTreeCoordinate parentCell(this->getX()>>1,this->getY()>>1,this->getZ()>>1);
@@ -373,7 +374,7 @@ public:
                             const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - this->getZ();
 
                             // Test if it is a direct neighbor
-                            if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                 // add to neighbors
                                 inNeighborsPosition[idxNeighbors] = ((( (xdiff+3) * 7) + (ydiff+3))) * 7 + zdiff + 3;
                                 inNeighbors[idxNeighbors++] = (mortonOther << 3) | idxCousin;
@@ -387,7 +388,7 @@ public:
         return idxNeighbors;
     }
 
-    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[189]) const{
+    int getInteractionNeighbors(const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], const int neighSeparation = 1) const{
         // Then take each child of the parent's neighbors if not in directNeighbors
         // Father coordinate
         const FTreeCoordinate parentCell(this->getX()>>1,this->getY()>>1,this->getZ()>>1);
@@ -418,7 +419,7 @@ public:
                             const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - this->getZ();
 
                             // Test if it is a direct neighbor
-                            if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                 // add to neighbors
                                 inNeighbors[idxNeighbors++] = (mortonOther << 3) | idxCousin;
                             }
diff --git a/Src/Core/FFmmAlgorithm.hpp b/Src/Core/FFmmAlgorithm.hpp
index 5526e19af..3b102ddc6 100644
--- a/Src/Core/FFmmAlgorithm.hpp
+++ b/Src/Core/FFmmAlgorithm.hpp
@@ -47,6 +47,7 @@ class FFmmAlgorithm :  public FAbstractAlgorithm, public FAlgorithmTimers {
 
     const int OctreeHeight;        ///< The height of the given tree.
 
+    const int leafLevelSeperationCriteria;
 public:
     /** Class constructor
      * 
@@ -56,8 +57,8 @@ public:
      *
      * \except An exception is thrown if one of the arguments is NULL.
      */
-    FFmmAlgorithm(OctreeClass* const inTree, KernelClass* const inKernels)
-        : tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()) {
+    FFmmAlgorithm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
+        : tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
         FAssertLF(kernels, "kernels cannot be null");
@@ -190,13 +191,16 @@ protected:
 
         const CellClass* neighbors[343];
 
+
         // for each levels
         for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
 
+            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
+
             // for each cells
             do{
-                const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel);
+                const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                 FLOG(computationCounter.tic());
                 if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, idxLevel);
                 FLOG(computationCounter.tac());
diff --git a/Src/Core/FFmmAlgorithmPeriodic.hpp b/Src/Core/FFmmAlgorithmPeriodic.hpp
index e82dfce6a..bdd5f8cca 100644
--- a/Src/Core/FFmmAlgorithmPeriodic.hpp
+++ b/Src/Core/FFmmAlgorithmPeriodic.hpp
@@ -51,6 +51,7 @@ class FFmmAlgorithmPeriodic : public FAbstractAlgorithm{
     const int nbLevelsAboveRoot;    //< The nb of level the user ask to go above the tree (>= -1)
     const int offsetRealTree;       //< nbLevelsAboveRoot GetFackLevel
 
+    const int leafLevelSeperationCriteria;
 
 public:
     /** The constructor need the octree and the kernels used for computation
@@ -60,9 +61,9 @@ public:
       * @param inUpperLevel this parameter defines the behavior of the periodicity refer to the main doc
       *
       */
-    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0)
+    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(nullptr), OctreeHeight(tree->getHeight()),
-          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3) {
+          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
         FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
@@ -254,9 +255,10 @@ protected:
         for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
             const int fackLevel = idxLevel + offsetRealTree;
+            const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
             // for each cells
             do{
-                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs);
+                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
                 FLOG(computationCounter.tic());
                 if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
                 FLOG(computationCounter.tac());
diff --git a/Src/Core/FFmmAlgorithmSectionTask.hpp b/Src/Core/FFmmAlgorithmSectionTask.hpp
index 0a797b499..0544d8481 100644
--- a/Src/Core/FFmmAlgorithmSectionTask.hpp
+++ b/Src/Core/FFmmAlgorithmSectionTask.hpp
@@ -49,15 +49,16 @@ class FFmmAlgorithmSectionTask : public FAbstractAlgorithm{
 
     const int OctreeHeight;
 
+    const int leafLevelSeperationCriteria;
 public:
     /** The constructor need the octree and the kernels used for computation
       * @param inTree the octree to work on
       * @param inKernels the kernels to call
       * An assert is launched if one of the arguments is null
       */
-    FFmmAlgorithmSectionTask(OctreeClass* const inTree, KernelClass* const inKernels)
+    FFmmAlgorithmSectionTask(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(0),
-          MaxThreads(omp_get_max_threads()), OctreeHeight(tree->getHeight())
+          MaxThreads(omp_get_max_threads()), OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria)
     {
 
         FAssertLF(tree, "tree cannot be null");
@@ -214,9 +215,10 @@ protected:
         // for each levels
         for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
+            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
             // for each cells
             do{
-                int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel);
+                const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                 if(counter){
                     #pragma omp task firstprivate(octreeIterator, neighbors, counter) shared(idxLevel)
                     {
diff --git a/Src/Core/FFmmAlgorithmTask.hpp b/Src/Core/FFmmAlgorithmTask.hpp
index 9301a76b9..84cd2c4c8 100644
--- a/Src/Core/FFmmAlgorithmTask.hpp
+++ b/Src/Core/FFmmAlgorithmTask.hpp
@@ -49,15 +49,16 @@ class FFmmAlgorithmTask : public FAbstractAlgorithm, public FAlgorithmTimers {
 
     const int OctreeHeight;
 
+    const int leafLevelSeperationCriteria;
 public:
     /** The constructor need the octree and the kernels used for computation
       * @param inTree the octree to work on
       * @param inKernels the kernels to call
       * An assert is launched if one of the arguments is null
       */
-    FFmmAlgorithmTask(OctreeClass* const inTree, KernelClass* const inKernels)
+    FFmmAlgorithmTask(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(nullptr),
-          MaxThreads(omp_get_max_threads()), OctreeHeight(tree->getHeight())
+          MaxThreads(omp_get_max_threads()), OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria)
     {
 
         FAssertLF(tree, "tree cannot be null");
@@ -228,9 +229,10 @@ protected:
                 // for each levels
                 for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
                     FLOG(FTic counterTimeLevel);
+                    const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
                     // for each cells
                     do{
-                        int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel);
+                        const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                         if(counter){
                             #pragma omp task firstprivate(octreeIterator, neighbors, counter) shared(idxLevel)
                             {
diff --git a/Src/Core/FFmmAlgorithmThread.hpp b/Src/Core/FFmmAlgorithmThread.hpp
index 529fc1488..5f850aaaf 100644
--- a/Src/Core/FFmmAlgorithmThread.hpp
+++ b/Src/Core/FFmmAlgorithmThread.hpp
@@ -62,6 +62,8 @@ class FFmmAlgorithmThread : public FAbstractAlgorithm, public FAlgorithmTimers{
 
     const bool staticSchedule;
 
+    const int leafLevelSeperationCriteria;
+
     template <class NumType>
     NumType getChunkSize(const NumType inSize) const {
         if(staticSchedule){
@@ -83,10 +85,10 @@ public:
      * \except An exception is thrown if one of the arguments is NULL.
      */
     FFmmAlgorithmThread(OctreeClass* const inTree, KernelClass* const inKernels,
-                        const bool inStaticSchedule = true)
+                        const bool inStaticSchedule = true, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(nullptr), iterArray(nullptr), leafsNumber(0),
           MaxThreads(omp_get_max_threads()), OctreeHeight(tree->getHeight()),
-          staticSchedule(inStaticSchedule) {
+          staticSchedule(inStaticSchedule), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
 
@@ -278,6 +280,7 @@ protected:
         // for each levels
         for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
+            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
             int numberOfCells = 0;
             // for each cells
             do{
@@ -297,7 +300,7 @@ protected:
 
                 #pragma omp for  schedule(dynamic, chunkSize) nowait
                 for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
-                    const int counter = tree->getInteractionNeighbors(neighbors,  iterArray[idxCell].getCurrentGlobalCoordinate(),idxLevel);
+                    const int counter = tree->getInteractionNeighbors(neighbors, iterArray[idxCell].getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                     if(counter) myThreadkernels->M2L( iterArray[idxCell].getCurrentCell() , neighbors, counter, idxLevel);
                 }
 
diff --git a/Src/Core/FFmmAlgorithmThreadProc.hpp b/Src/Core/FFmmAlgorithmThreadProc.hpp
index c96c035f2..2115736ac 100644
--- a/Src/Core/FFmmAlgorithmThreadProc.hpp
+++ b/Src/Core/FFmmAlgorithmThreadProc.hpp
@@ -81,6 +81,7 @@ class FFmmAlgorithmThreadProc : public FAbstractAlgorithm, public FAlgorithmTime
 
     const int OctreeHeight;            //<Height of the tree
 
+    const int leafLevelSeperationCriteria;
 
     /** An interval is the morton index interval
      * that a proc use (it holds data in this interval)
@@ -134,10 +135,12 @@ public:
      * @param inKernels the kernels to call
      * An assert is launched if one of the arguments is null
      */
-    FFmmAlgorithmThreadProc(const FMpi::FComm& inComm, OctreeClass* const inTree, KernelClass* const inKernels)
+    FFmmAlgorithmThreadProc(const FMpi::FComm& inComm, OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(nullptr), comm(inComm), iterArray(nullptr),iterArrayComm(nullptr),numberOfLeafs(0),
           MaxThreads(omp_get_max_threads()), nbProcess(inComm.processCount()), idProcess(inComm.processId()),
-          OctreeHeight(tree->getHeight()),intervals(new Interval[inComm.processCount()]),
+          OctreeHeight(tree->getHeight()),
+          leafLevelSeperationCriteria(inLeafLevelSeperationCriteria),
+          intervals(new Interval[inComm.processCount()]),
           workingIntervalsPerLevel(new Interval[inComm.processCount() * tree->getHeight()])
     {
         FAssertLF(tree, "tree cannot be null");
@@ -622,6 +625,9 @@ protected:
                     typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
                     // for each levels
                     for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
+
+                        const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
+
                         if(!procHasWorkAtLevel(idxLevel, idProcess)){
                             avoidGotoLeftIterator.moveDown();
                             octreeIterator = avoidGotoLeftIterator;
@@ -646,10 +652,10 @@ protected:
 
                         // Which cell potentialy needs other data and in the same time
                         // are potentialy needed by other
-                        MortonIndex neighborsIndexes[189];
+                        MortonIndex neighborsIndexes[/*189+26+1*/216];
                         for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                             // Find the M2L neigbors of a cell
-                            const int counter = iterArrayLocal[idxCell].getCurrentGlobalCoordinate().getInteractionNeighbors(idxLevel, neighborsIndexes);
+                            const int counter = iterArrayLocal[idxCell].getCurrentGlobalCoordinate().getInteractionNeighbors(idxLevel, neighborsIndexes, separationCriteria);
 
                             memset(alreadySent, false, sizeof(bool) * nbProcess);
                             bool needOther = false;
@@ -778,6 +784,8 @@ protected:
                 // Now we can compute all the data
                 // for each levels
                 for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
+                    const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
+
                     if(!procHasWorkAtLevel(idxLevel, idProcess)){
                         avoidGotoLeftIterator.moveDown();
                         octreeIterator = avoidGotoLeftIterator;
@@ -807,7 +815,7 @@ protected:
 
                                 const int nbCellToCompute = FMath::Min(chunckSize, numberOfCells-idxCell);
                                 for(int idxCellToCompute = idxCell ; idxCellToCompute < idxCell+nbCellToCompute ; ++idxCellToCompute){
-                                    const int counter = tree->getInteractionNeighbors(neighbors,  iterArray[idxCellToCompute].getCurrentGlobalCoordinate(), idxLevel);
+                                    const int counter = tree->getInteractionNeighbors(neighbors,  iterArray[idxCellToCompute].getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                                     if(counter) myThreadkernels->M2L( iterArray[idxCellToCompute].getCurrentCell() , neighbors, counter, idxLevel);
                                 }
                             }
@@ -843,6 +851,8 @@ protected:
             // compute the second time
             // for each levels
             for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
+                const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
+
                 if(!procHasWorkAtLevel(idxLevel, idProcess)){
                     avoidGotoLeftIterator.moveDown();
                     octreeIterator = avoidGotoLeftIterator;
@@ -900,15 +910,15 @@ protected:
                 #pragma omp parallel
                 {
                     KernelClass * const myThreadkernels = kernels[omp_get_thread_num()];
-                    MortonIndex neighborsIndex[189];
-                    int neighborsPosition[189];
+                    MortonIndex neighborsIndex[/*189+26+1*/216];
+                    int neighborsPosition[/*189+26+1*/216];
                     const CellClass* neighbors[343];
 
                     #pragma omp for schedule(static) nowait
                     for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                         // compute indexes
                         memset(neighbors, 0, 343 * sizeof(CellClass*));
-                        const int counterNeighbors = iterArray[idxCell].getCurrentGlobalCoordinate().getInteractionNeighbors(idxLevel, neighborsIndex, neighborsPosition);
+                        const int counterNeighbors = iterArray[idxCell].getCurrentGlobalCoordinate().getInteractionNeighbors(idxLevel, neighborsIndex, neighborsPosition, separationCriteria);
 
                         int counter = 0;
                         // does we receive this index from someone?
diff --git a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
index 54605f98e..a68cd09eb 100644
--- a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
+++ b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp
@@ -83,6 +83,8 @@ class FFmmAlgorithmThreadProcPeriodic : public FAbstractAlgorithm {
 
     const int OctreeHeight;
 
+    const int leafLevelSeperationCriteria;
+
 public:
     struct Interval{
         MortonIndex leftIndex;
@@ -139,12 +141,13 @@ public:
      * An assert is launched if one of the arguments is null
      */
     FFmmAlgorithmThreadProcPeriodic(const FMpi::FComm& inComm, OctreeClass* const inTree,
-                                    const int inUpperLevel = 2)
+                                    const int inUpperLevel = 2, const int inLeafLevelSeperationCriteria = 1)
         : tree(inTree) , kernels(nullptr), comm(inComm), nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3),
           numberOfLeafs(0),
           MaxThreads(omp_get_max_threads()), nbProcess(inComm.processCount()), idProcess(inComm.processId()),
           OctreeHeight(tree->getHeight()),intervals(new Interval[inComm.processCount()]),
-          workingIntervalsPerLevel(new Interval[inComm.processCount() * tree->getHeight()]) {
+          workingIntervalsPerLevel(new Interval[inComm.processCount() * tree->getHeight()]),
+          leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
         FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
@@ -752,6 +755,9 @@ protected:
                     typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
                     // for each levels
                     for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
+
+                        const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
+
                         if(!procHasWorkAtLevel(idxLevel, idProcess)){
                             avoidGotoLeftIterator.moveDown();
                             octreeIterator = avoidGotoLeftIterator;
@@ -776,13 +782,13 @@ protected:
 
                         // Which cell potentialy needs other data and in the same time
                         // are potentialy needed by other
-                        int neighborsPosition[189];
-                        MortonIndex neighborsIndexes[189];
+                        int neighborsPosition[/*189+26+1*/216];
+                        MortonIndex neighborsIndexes[/*189+26+1*/216];
                         for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                             // Find the M2L neigbors of a cell
                             const int counter = getPeriodicInteractionNeighbors(iterArray[idxCell].getCurrentGlobalCoordinate(),
                                                                                idxLevel,
-                                                                               neighborsIndexes, neighborsPosition, AllDirs);
+                                                                               neighborsIndexes, neighborsPosition, AllDirs, leafLevelSeperationCriteria);
 
                             memset(alreadySent, false, sizeof(bool) * nbProcess);
                             bool needOther = false;
@@ -907,6 +913,9 @@ protected:
                 // for each levels
                 for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
                     const int fackLevel = idxLevel + offsetRealTree;
+
+                    const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
+
                     if(!procHasWorkAtLevel(idxLevel, idProcess)){
                         avoidGotoLeftIterator.moveDown();
                         octreeIterator = avoidGotoLeftIterator;
@@ -936,10 +945,9 @@ protected:
 
                                 const int nbCellToCompute = FMath::Min(chunckSize, numberOfCells-idxCell);
                                 for(int idxCellToCompute = idxCell ; idxCellToCompute < idxCell+nbCellToCompute ; ++idxCellToCompute){
-                                    const int counter = tree->
-                                            getPeriodicInteractionNeighbors(neighbors,
+                                    const int counter = tree->getPeriodicInteractionNeighbors(neighbors,
                                                                             iterArray[idxCellToCompute].getCurrentGlobalCoordinate(),
-                                                                            idxLevel, AllDirs);
+                                                                            idxLevel, AllDirs, separationCriteria);
 
                                     if(counter)
                                         myThreadkernels->M2L( iterArray[idxCellToCompute].getCurrentCell() , neighbors, counter, fackLevel);
@@ -973,6 +981,9 @@ protected:
             // for each levels
             for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
                 const int fackLevel = idxLevel + offsetRealTree;
+
+                const int separationCriteria = (fackLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
+
                 if(!procHasWorkAtLevel(idxLevel, idProcess)){
                     avoidGotoLeftIterator.moveDown();
                     octreeIterator = avoidGotoLeftIterator;
@@ -1030,15 +1041,15 @@ protected:
                 #pragma omp parallel
                 {
                     KernelClass * const myThreadkernels = kernels[omp_get_thread_num()];
-                    MortonIndex neighborsIndex[189];
-                    int neighborsPosition[189];
+                    MortonIndex neighborsIndex[/*189+26+1*/216];
+                    int neighborsPosition[/*189+26+1*/216];
                     const CellClass* neighbors[343];
 
                     #pragma omp for schedule(static) nowait
                     for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                         // compute indexes
                         memset(neighbors, 0, 343 * sizeof(CellClass*));
-                        const int counterNeighbors = getPeriodicInteractionNeighbors(iterArray[idxCell].getCurrentGlobalCoordinate(), idxLevel, neighborsIndex, neighborsPosition, AllDirs);
+                        const int counterNeighbors = getPeriodicInteractionNeighbors(iterArray[idxCell].getCurrentGlobalCoordinate(), idxLevel, neighborsIndex, neighborsPosition, AllDirs, separationCriteria);
                         //const int counterNeighbors = iterArray[idxCell].getCurrentGlobalCoordinate().getInteractionNeighbors(idxLevel, neighborsIndex, neighborsPosition);
 
                         int counter = 0;
@@ -1758,7 +1769,8 @@ protected:
     }
 
 
-    int getPeriodicInteractionNeighbors(const FTreeCoordinate& workingCell,const int inLevel, MortonIndex inNeighbors[189], int inNeighborsPosition[189], const int inDirection) const{
+    int getPeriodicInteractionNeighbors(const FTreeCoordinate& workingCell,const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], int inNeighborsPosition[/*189+26+1*/216],
+                const int inDirection, const int neighSeparation) const{
 
         // Then take each child of the parent's neighbors if not in directNeighbors
         // Father coordinate
@@ -1821,7 +1833,7 @@ protected:
                             const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1))         - workingCell.getZ();
 
                             // Test if it is a direct neighbor
-                            if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                 // add to neighbors
                                 inNeighborsPosition[idxNeighbors] = (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3;
                                 inNeighbors[idxNeighbors++] = (mortonOtherParent << 3) | idxCousin;
@@ -1835,7 +1847,8 @@ protected:
         return idxNeighbors;
     }
 
-    int getInteractionNeighbors(const FTreeCoordinate& workingCell,const int inLevel, MortonIndex inNeighbors[189], int inNeighborsPosition[189]) const{
+    int getInteractionNeighbors(const FTreeCoordinate& workingCell,const int inLevel, MortonIndex inNeighbors[/*189+26+1*/216], int inNeighborsPosition[/*189+26+1*/216],
+                const int neighSeparation) const{
 
         // Then take each child of the parent's neighbors if not in directNeighbors
         // Father coordinate
@@ -1860,7 +1873,7 @@ protected:
                             const int zdiff  = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
 
                             // Test if it is a direct neighbor
-                            if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+                            if(FMath::Abs(xdiff) > neighSeparation || FMath::Abs(ydiff) > neighSeparation || FMath::Abs(zdiff) > neighSeparation){
                                 // add to neighbors
                                 inNeighborsPosition[idxNeighbors] = ((( (xdiff+3) * 7) + (ydiff+3))) * 7 + zdiff + 3;
                                 inNeighbors[idxNeighbors++] = (mortonOtherParent << 3) | idxCousin;
diff --git a/Src/Core/FFmmAlgorithmThreadTsm.hpp b/Src/Core/FFmmAlgorithmThreadTsm.hpp
index 86434c674..6bb6a268a 100644
--- a/Src/Core/FFmmAlgorithmThreadTsm.hpp
+++ b/Src/Core/FFmmAlgorithmThreadTsm.hpp
@@ -55,15 +55,17 @@ class FFmmAlgorithmThreadTsm : public FAbstractAlgorithm, public FAlgorithmTimer
 
     const int OctreeHeight;
 
+    const int leafLevelSeperationCriteria;
+
 public:
     /** The constructor need the octree and the kernels used for computation
       * @param inTree the octree to work on
       * @param inKernels the kernels to call
       * An assert is launched if one of the arguments is null
       */
-    FFmmAlgorithmThreadTsm(OctreeClass* const inTree, KernelClass* const inKernels)
+    FFmmAlgorithmThreadTsm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
                       : tree(inTree) , kernels(nullptr), iterArray(nullptr),
-                      MaxThreads(omp_get_max_threads()) , OctreeHeight(tree->getHeight()) {
+                      MaxThreads(omp_get_max_threads()) , OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
 
@@ -248,6 +250,7 @@ protected:
             // for each levels
             for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
                 FLOG(FTic counterTimeLevel);
+                const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
 
                 int numberOfCells = 0;
                 // for each cells
@@ -270,7 +273,7 @@ protected:
                     for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                         CellClass* const currentCell = iterArray[idxCell].getCurrentCell();
                         if(currentCell->hasTargetsChild()){
-                            const int counter = tree->getInteractionNeighbors(neighbors, iterArray[idxCell].getCurrentGlobalCoordinate(),idxLevel);
+                            const int counter = tree->getInteractionNeighbors(neighbors, iterArray[idxCell].getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
                             if( counter ){
                                 int counterWithSrc = 0;
                                 for(int idxRealNeighbors = 0 ; idxRealNeighbors < 343 ; ++idxRealNeighbors ){
diff --git a/Src/Core/FFmmAlgorithmTsm.hpp b/Src/Core/FFmmAlgorithmTsm.hpp
index 2937a880b..80f238d68 100644
--- a/Src/Core/FFmmAlgorithmTsm.hpp
+++ b/Src/Core/FFmmAlgorithmTsm.hpp
@@ -46,6 +46,8 @@ class FFmmAlgorithmTsm : public FAbstractAlgorithm{
 
     const int OctreeHeight;
 
+    const int leafLevelSeperationCriteria;
+
     FLOG(FTic counterTime);                                               //< In case of debug: to count the elapsed time
     FLOG(FTic computationCounter);                                        //< In case of debug: to  count computation time
 
@@ -55,8 +57,8 @@ public:
       * @param inKernels the kernels to call
       * An assert is launched if one of the arguments is null
       */
-    FFmmAlgorithmTsm(OctreeClass* const inTree, KernelClass* const inKernels)
-        : tree(inTree) , kernels(inKernels) , OctreeHeight(tree->getHeight()){
+    FFmmAlgorithmTsm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
+        : tree(inTree) , kernels(inKernels) , OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria){
 
         FAssertLF(tree, "tree cannot be null");
         FAssertLF(kernels, "kernels cannot be null");
@@ -198,13 +200,14 @@ protected:
         // for each levels
         for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
+            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
             // for each cells
             do{
                 FLOG(computationCounter.tic());
                 CellClass* const currentCell = octreeIterator.getCurrentCell();
 
                 if(currentCell->hasTargetsChild()){
-                    const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),idxLevel);
+                    const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),idxLevel, separationCriteria);
                     if( counter ){
                         int counterWithSrc = 0;
                         for(int idxRealNeighbors = 0 ; idxRealNeighbors < 343 ; ++idxRealNeighbors ){
-- 
GitLab