FFmmAlgorithmPeriodic.hpp 58.8 KB
Newer Older
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15
// ===================================================================================
16 17
#ifndef FFMMALGORITHMPERIODIC_HPP
#define FFMMALGORITHMPERIODIC_HPP
18

19 20

#include "../Utils/FGlobal.hpp"
21
#include "../Utils/FGlobalPeriodic.hpp"
22 23 24 25
#include "../Utils/FAssertable.hpp"
#include "../Utils/FDebug.hpp"
#include "../Utils/FTrace.hpp"
#include "../Utils/FTic.hpp"
26
#include "../Utils/FMemUtils.hpp"
27 28 29 30

#include "../Containers/FOctree.hpp"
#include "../Containers/FVector.hpp"

31
#include "FCoreCommon.hpp"
32 33 34 35 36 37 38

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFmmAlgorithmPeriodic
* @brief
* Please read the license
*
39
* This class is a basic FMM algorithm with periodic behavior
40 41 42 43
* It just iterates on a tree and call the kernels with good arguments.
*
* Of course this class does not deallocate pointer given in arguements.
*/
44
template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
45 46
class FFmmAlgorithmPeriodic : protected FAssertable{

47 48
    OctreeClass* const tree;        //< The octree to work on
    KernelClass* kernels;           //< The kernels
49

50 51 52 53 54
    const int OctreeHeight;         //< The heigh of the octree (real heigh)
    const int nbLevelsAboveRoot;    //< The nb of level the user ask to go above the tree (>= -1)
    const int offsetRealTree;       //< nbLevelsAboveRoot GetFackLevel
    const int periodicDirections;

55 56 57 58 59 60

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
61 62
      * @param inUpperLevel this parameter defins the behavior of the periodicity refer to the main doc
      *
63
      */
64 65
    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inPeriodicDirections = AllDirs)
        : tree(inTree) , kernels(0), OctreeHeight(tree->getHeight()),
66
          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3),
67
          periodicDirections(inPeriodicDirections) {
68 69

        fassert(tree, "tree cannot be null", __LINE__, __FILE__);
70
        fassert(-1 <= inUpperLevel, "inUpperLevel cannot be < -1", __LINE__, __FILE__);
71 72 73 74 75 76 77 78

        FDEBUG(FDebug::Controller << "FFmmAlgorithmPeriodic\n");
    }

    /** Default destructor */
    virtual ~FFmmAlgorithmPeriodic(){
    }

79 80 81 82
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

83 84 85 86
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
87
    void execute(const unsigned operationsToProceed = FFmmNearAndFarFields){
88
        fassert(kernels, "kernels cannot be null", __LINE__, __FILE__);
89
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );       
90

91
        if(operationsToProceed & FFmmP2M) bottomPass();
92

93
        if(operationsToProceed & FFmmM2M) upwardPass();
94

95 96 97 98 99
        if(operationsToProceed & FFmmM2L){
            transferPass();
            // before downward pass we have to perform the periodicity
            processPeriodicLevels();
        }
100

101
        if(operationsToProceed & FFmmL2L) downardPass();
102

103
        if(operationsToProceed & FFmmP2P || operationsToProceed & FFmmL2P) directPass();
104 105
    }

106

107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
        FDEBUG( FDebug::Controller.write("\tStart Bottom Pass\n").write(FDebug::Flush) );
        FDEBUG(FTic counterTime);
        FDEBUG(FTic computationCounter);

        typename OctreeClass::Iterator octreeIterator(tree);

        // Iterate on leafs
        octreeIterator.gotoBottomLeft();
        do{
            // We need the current cell that represent the leaf
            // and the list of particles
            FDEBUG(computationCounter.tic());
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
            FDEBUG(computationCounter.tac());
        } while(octreeIterator.moveRight());

        FDEBUG( FDebug::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
    }

    /////////////////////////////////////////////////////////////////////////////
    // Upward
    /////////////////////////////////////////////////////////////////////////////

    /** M2M */
    void upwardPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
        FDEBUG( FDebug::Controller.write("\tStart Upward Pass\n").write(FDebug::Flush); );
        FDEBUG(FTic counterTime);
        FDEBUG(FTic computationCounter);

        // Start from leal level - 1
        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        octreeIterator.moveUp();

        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);

        // for each levels
        for(int idxLevel = OctreeHeight - 2 ; idxLevel > 0 ; --idxLevel ){
154
            FDEBUG(FTic counterTimeLevel);
155
            const int fackLevel = idxLevel + offsetRealTree;
156 157 158 159 160
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
                FDEBUG(computationCounter.tic());
161
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
162 163 164 165 166
                FDEBUG(computationCounter.tac());
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
167
            FDEBUG( FDebug::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
168 169 170 171 172 173 174 175
        }


        FDEBUG( FDebug::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
    }

    /////////////////////////////////////////////////////////////////////////////
176
    // Transfer
177 178 179
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
180
    void transferPass(){
181 182
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );

183 184 185 186 187 188 189
        FDEBUG( FDebug::Controller.write("\tStart Downward Pass (M2L)\n").write(FDebug::Flush); );
        FDEBUG(FTic counterTime);
        FDEBUG(FTic computationCounter);

        typename OctreeClass::Iterator octreeIterator(tree);
        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);

190
        const CellClass* neighbors[343];
191 192 193

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
194
            FDEBUG(FTic counterTimeLevel);
195
            const int fackLevel = idxLevel + offsetRealTree;
196 197
            // for each cells
            do{
198
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, periodicDirections);
199
                FDEBUG(computationCounter.tic());
200
                if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
201 202 203 204
                FDEBUG(computationCounter.tac());
            } while(octreeIterator.moveRight());
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
205 206

            FDEBUG(computationCounter.tic());
207
            kernels->finishedLevelM2L(fackLevel);
208
            FDEBUG(computationCounter.tac());
209
            FDEBUG( FDebug::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
210
        }
211 212 213
        FDEBUG( FDebug::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
    }
214

215 216 217 218 219 220 221 222 223 224 225 226 227
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
        FDEBUG( FDebug::Controller.write("\tStart Downward Pass (L2L)\n").write(FDebug::Flush); );
        FDEBUG(FTic counterTime);
        FDEBUG(FTic computationCounter );

        typename OctreeClass::Iterator octreeIterator(tree);
        typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
228

229 230 231
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
232
            FDEBUG(FTic counterTimeLevel);
233 234
            const int fackLevel = idxLevel + offsetRealTree;

235 236 237
            // for each cells
            do{
                FDEBUG(computationCounter.tic());
238
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
239 240 241 242 243
                FDEBUG(computationCounter.tac());
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
244
            FDEBUG( FDebug::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
245 246
        }

247 248 249
        FDEBUG( FDebug::Controller << "\tFinished (@Downward Pass (L2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );

250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265

    }

    /////////////////////////////////////////////////////////////////////////////
    // Direct
    /////////////////////////////////////////////////////////////////////////////

    /** P2P */
    void directPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
        FDEBUG( FDebug::Controller.write("\tStart Direct Pass\n").write(FDebug::Flush); );
        FDEBUG(FTic counterTime);
        FDEBUG(FTic computationCounterL2P);
        FDEBUG(FTic computationCounterP2P);

        const int heightMinusOne = OctreeHeight - 1;
266
        const FReal boxWidth = tree->getBoxWidth();
267 268 269 270

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
berenger-bramas's avatar
berenger-bramas committed
271
        ContainerClass* neighbors[27];
272
        FTreeCoordinate offsets[27];
273
        bool hasPeriodicLeaves;
274 275 276 277 278
        // for each leafs
        do{
            FDEBUG(computationCounterL2P.tic());
            kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
            FDEBUG(computationCounterL2P.tac());
279

280
            // need the current particles and neighbors particles
281
            const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
282 283 284 285 286 287
            const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, periodicDirections);
            int periodicNeighborsCounter = 0;

            if(hasPeriodicLeaves){
                ContainerClass* periodicNeighbors[27];
                memset(periodicNeighbors, 0, 27 * sizeof(ContainerClass*));
288 289

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
290 291 292 293 294
                    if( neighbors[idxNeig] && !offsets[idxNeig].equals(0,0,0) ){
                        // Put periodic neighbors into other array
                        periodicNeighbors[idxNeig] = neighbors[idxNeig];
                        neighbors[idxNeig] = 0;
                        ++periodicNeighborsCounter;
295 296 297 298 299 300 301 302 303

                        FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
                        FReal*const positionsY = periodicNeighbors[idxNeig]->getWPositions()[1];
                        FReal*const positionsZ = periodicNeighbors[idxNeig]->getWPositions()[2];

                        for(int idxPart = 0; idxPart < periodicNeighbors[idxNeig]->getNbParticles() ; ++idxPart){
                            positionsX[idxPart] += boxWidth * FReal(offsets[idxNeig].getX());
                            positionsY[idxPart] += boxWidth * FReal(offsets[idxNeig].getY());
                            positionsZ[idxPart] += boxWidth * FReal(offsets[idxNeig].getZ());
304 305 306 307
                        }
                    }
                }

308 309 310 311
                FDEBUG(computationCounterP2P.tic());
                kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                             octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborsCounter);
                FDEBUG(computationCounterP2P.tac());
312 313

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
314 315 316 317 318 319 320 321 322
                    if( periodicNeighbors[idxNeig] ){                        
                        FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
                        FReal*const positionsY = periodicNeighbors[idxNeig]->getWPositions()[1];
                        FReal*const positionsZ = periodicNeighbors[idxNeig]->getWPositions()[2];

                        for(int idxPart = 0; idxPart < periodicNeighbors[idxNeig]->getNbParticles() ; ++idxPart){
                            positionsX[idxPart] -= boxWidth * FReal(offsets[idxNeig].getX());
                            positionsY[idxPart] -= boxWidth * FReal(offsets[idxNeig].getY());
                            positionsZ[idxPart] -= boxWidth * FReal(offsets[idxNeig].getZ());
323 324 325 326 327
                        }
                    }
                }
            }

328 329 330 331 332 333
            FDEBUG(computationCounterP2P.tic());
            kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                         octreeIterator.getCurrentListSrc(), neighbors, counter - periodicNeighborsCounter);
            FDEBUG(computationCounterP2P.tac());


334 335 336 337 338 339 340 341 342 343 344 345 346
        } while(octreeIterator.moveRight());


        FDEBUG( FDebug::Controller << "\tFinished (@Direct Pass (L2P + P2P) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation L2P : " << computationCounterL2P.cumulated() << " s\n" );
        FDEBUG( FDebug::Controller << "\t\t Computation P2P : " << computationCounterP2P.cumulated() << " s\n" );

    }

    /////////////////////////////////////////////////////////////////////////////
    // Periodic levels = levels <= 0
    /////////////////////////////////////////////////////////////////////////////

347

348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
    /** This function process several M2M from level nbLevelsAboveRoot to level 0
      * and give the final result
      * @param result the cell at the last M2M
      * @param root the starting cell
      * @param startX the beginning of the index in x [0;endX]
      * @param endX the end of the index in x [startX;1]
      * @param startY the beginning of the index in y [0;endY]
      * @param endY the end of the index in y [startY;1]
      * @param startZ the beginning of the index in z [0;endZ]
      * @param endZ the end of the index in z [startZ;1]
      */
    void processTopM2MInIntervals( CellClass*const result, const CellClass& root, const int startX,
                              const int endX, const int startY, const int endY, const int startZ,
                              const int endZ){
        // allocate array
        CellClass*const cellsAtLevel = new CellClass[nbLevelsAboveRoot+2];
        // process by using other function
        processM2MInIntervals(cellsAtLevel,root,startX,endX,startY,endY,startZ,endZ);
        // copy result
        *result = cellsAtLevel[0];
        delete[] cellsAtLevel;
    }

    /** This function process several M2M from level nbLevelsAboveRoot to level 0
      * @param cellsAtLevel the intermediate results
      * @param root the starting cell
      * @param startX the beginning of the index in x [0;endX]
      * @param endX the end of the index in x [startX;1]
      * @param startY the beginning of the index in y [0;endY]
      * @param endY the end of the index in y [startY;1]
      * @param startZ the beginning of the index in z [0;endZ]
      * @param endZ the end of the index in z [startZ;1]
      */
    void  processM2MInIntervals( CellClass cellsAtLevel[], const CellClass& root, const int startX,
                              const int endX, const int startY, const int endY, const int startZ,
                              const int endZ){
        // start from the initial cell
        cellsAtLevel[nbLevelsAboveRoot+1] = root;
        // to create virtual children
        CellClass* virtualChild[8];
        // for all levels
        for(int idxLevel = nbLevelsAboveRoot ; idxLevel >= 0  ; --idxLevel){
            // reset children
            memset(virtualChild, 0, sizeof(CellClass*)*8);
            // fill the vector with previous result
            for(int idxX = startX ; idxX <= endX ; ++idxX){
                for(int idxY = startY ; idxY <= endY ; ++idxY){
                    for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
                        virtualChild[childIndex(idxX,idxY,idxZ)] = &cellsAtLevel[idxLevel+1];
                    }
                }
            }
            // compute the M2M
            kernels->M2M( &cellsAtLevel[idxLevel], virtualChild, idxLevel + 2);
        }
    }

    /** Fill an interactions neighbors with some intervals
      * @param neighbors the vector to fill
      * @param source the source cell to fill the vector
      * @param startX the beginning of the index in x [-3;0]
      * @param endX the end of the index in x  [0;3]
      * @param startY the beginning of the index in y [-3;0]
      * @param endY the end of the index in y [0;3]
      * @param startZ the beginning of the index in z [-3;0]
      * @param endZ the end of the index in z [0;3]
      * @return the number of position filled
      */
    int  fillM2LVectorFromIntervals(const CellClass* neighbors[343], const CellClass& source,
                     const int startX, const int endX, const int startY, const int endY,
                     const int startZ, const int endZ){
        int counter = 0;
        // for all x in interval
        for(int idxX = startX ; idxX <= endX ; ++idxX){
            // for all y in interval
            for(int idxY = startY ; idxY <= endY ; ++idxY){
                // for all z in interval
                for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
                    // do not fill close neigbors
427
                    if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1 ){
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457
                        neighbors[neighIndex(idxX,idxY,idxZ)] = &source;
                        ++counter;
                    }
                }
            }
        }
        // return the number of position filled
        return counter;
    }

    /** Get the index of a child (for the M2M and the L2L)
      * @param x the x position in the children  (from 0 to +1)
      * @param y the y position in the children  (from 0 to +1)
      * @param z the z position in the children  (from 0 to +1)
      * @return the index (from 0 to 7)
      */
    int childIndex(const int x, const int y, const int z) const {
        return (x<<2) | (y<<1) | z;
    }

    /** Get the index of a interaction neighbors (for M2L)
      * @param x the x position in the interactions (from -3 to +3)
      * @param y the y position in the interactions (from -3 to +3)
      * @param z the z position in the interactions (from -3 to +3)
      * @return the index (from 0 to 342)
      */
    int neighIndex(const int x, const int y, const int z) const {
        return (((x+3)*7) + (y+3))*7 + (z + 3);
    }

458 459
    /** To know how many times the box is repeated in each direction
      * -x +x -y +y -z +z
460 461 462
      * If the periodicy is not in all direction this number is unchanged
      * because it contains the theorical periodic box width for the
      * nbLevelsAboveRoot choosen
463
      */
464 465
    long long int theoricalRepetition() const {
        return nbLevelsAboveRoot == -1 ? 3 : 3LL * (1LL<<(nbLevelsAboveRoot+1)) + 1LL;
466 467
    }

468 469 470 471 472 473 474
    /** To know the number of box repeated in each direction
      * @param min the number of repetition for -x,-y,-z
      * @param max the number of repetition for x,y,z
      * The mins value are contains between [-(theoricalRepetition-1 / 2); 0]
      * The maxs value are contains between [0;(theoricalRepetition-1 / 2)]
      */
    void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
475
        const int halfRepeated = int((theoricalRepetition()-1)/2);
476 477
        min->setPosition(-ifDir(DirMinusX,halfRepeated,0),-ifDir(DirMinusY,halfRepeated,0),
                         -ifDir(DirMinusZ,halfRepeated,0));
478 479 480 481
        max->setPosition(ifDir(DirPlusX,halfRepeated,0),ifDir(DirPlusY,halfRepeated,0),
                         ifDir(DirPlusZ,halfRepeated,0));
    }

482 483 484 485
    /** To get the number of repetition in all direction (x,y,z)
      * @return the number of box in all directions
      * Each value is between [1;theoricalRepetition]
      */
486
    FTreeCoordinate repetitions() const {
487
        const int halfRepeated = int((theoricalRepetition()-1)/2);
488 489 490 491 492 493 494 495 496 497 498
        return FTreeCoordinate(ifDir(DirMinusX,halfRepeated,0) + ifDir(DirPlusX,halfRepeated,0) + 1,
                               ifDir(DirMinusY,halfRepeated,0) + ifDir(DirPlusY,halfRepeated,0) + 1,
                               ifDir(DirMinusZ,halfRepeated,0) + ifDir(DirPlusZ,halfRepeated,0) + 1);
    }

    /** This function has to be used to init the kernel with correct args
      * it return the box seen from a kernel point of view from the periodicity the user ask for
      * this is computed using the originalBoxWidth given in parameter
      * @param originalBoxWidth the real system size
      * @return the size the kernel should use
      */
499 500
    FReal extendedBoxWidth() const {
        return tree->getBoxWidth() * FReal(1<<offsetRealTree);
501 502
    }

503 504 505 506 507 508 509 510 511
    /** This function has to be used to init the kernel with correct args
      * it return the box cneter seen from a kernel point of view from the periodicity the user ask for
      * this is computed using the originalBoxWidth and originalBoxCenter given in parameter
      * @param originalBoxCenter the real system center
      * @param originalBoxWidth the real system size
      * @return the center the kernel should use
      */
    FPoint extendedBoxCenter() const {
        const FReal originalBoxWidth = tree->getBoxWidth();
512
        const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
513
        const FPoint originalBoxCenter = tree->getBoxCenter();
514 515 516 517 518

        const FReal offset = extendedBoxWidth()/2;
        return FPoint( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
519 520 521 522 523 524 525 526
    }

    /** This function has to be used to init the kernel with correct args
      * it return the tree heigh seen from a kernel point of view from the periodicity the user ask for
      * this is computed using the originalTreeHeight given in parameter
      * @param originalTreeHeight the real tree heigh
      * @return the heigh the kernel should use
      */
527
    int extendedTreeHeight() const {
528 529 530 531
        // The real height
        return OctreeHeight + offsetRealTree;
    }

532 533 534 535
    /** To know if a direction is used
      * @param testDir a direction to test
      * @return true if the direction is used else false
      */
536
    bool usePerDir(const int testDir) const{
537
        return TestPeriodicCondition(periodicDirections , PeriodicCondition(testDir));
538 539
    }

540 541 542 543 544 545
    /** To enable quick test of the direction
      * @param testDir the direction to test
      * @param correctValue the value to return if direction is use
      * @param wrongValue the value to return if direction is not use
      * @return correctValue if testDir is used, else wrongValue
      */
546
    template <class T>
547
    const T& ifDir(const PeriodicCondition testDir, const T& correctValue, const T& wrongValue) const {
548 549 550
        return (periodicDirections & testDir ? correctValue : wrongValue);
    }

551 552 553 554 555 556 557 558 559 560 561
    /** Periodicity Core
      * This function is split in several part:
      * 1 - special case managment
      * There is nothing to do if nbLevelsAboveRoot == -1 and only
      * a M2L if nbLevelsAboveRoot == 0
      * 2 - if nbLevelsAboveRoot > 0
      * First we compute M2M and special M2M if needed for the border
      * Then the M2L by taking into account the periodicity directions
      * Then the border by using the precomputed M2M
      * Finally the L2L
      */
562
    void processPeriodicLevels(){
563 564 565
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
        FDEBUG( FDebug::Controller.write("\tStart Periodic Pass\n").write(FDebug::Flush); );
        FDEBUG(FTic counterTime);
566
        /////////////////////////////////////////////////////
567 568
        // If nb level == -1 nothing to do
        if( nbLevelsAboveRoot == -1 ){
569
            FDEBUG( FDebug::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
570 571
            return;
        }
572
        /////////////////////////////////////////////////////
573 574 575 576 577 578
        // if nb level == 0 only M2L at real root level
        if( nbLevelsAboveRoot == 0 ){
            CellClass rootUp;
            // compute the root
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoLeft();
579
            kernels->M2M( &rootUp, octreeIterator.getCurrentBox(), 3);
580 581 582 583 584

            // build fack M2L vector from -3/+3 x/y/z
            const CellClass* neighbors[343];
            memset(neighbors, 0, sizeof(CellClass*) * 343);
            int counter = 0;
585 586 587
            for(int idxX = ifDir(DirMinusX,-3,0) ; idxX <= ifDir(DirPlusX,3,0) ; ++idxX){
                for(int idxY = ifDir(DirMinusY,-3,0) ; idxY <= ifDir(DirPlusY,3,0) ; ++idxY){
                    for(int idxZ = ifDir(DirMinusZ,-3,0) ; idxZ <= ifDir(DirPlusZ,3,0) ; ++idxZ){
588
                        if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
589
                            neighbors[neighIndex(idxX,idxY,idxZ)] = &rootUp;
590 591 592 593 594 595 596
                            ++counter;
                        }
                    }
                }
            }
            // compute M2L
            CellClass rootDown;
597
            kernels->M2L( &rootDown , neighbors, counter, 3);
598 599

            // put result in level 1
600
            kernels->L2L( &rootDown, octreeIterator.getCurrentBox(), 3);
601

602
            FDEBUG( FDebug::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
603 604
            return;
        }
605
        /////////////////////////////////////////////////////
606 607
        // in other situation, we have to compute M2L from 0 to nbLevelsAboveRoot
        // but also at nbLevelsAboveRoot +1 for the rest
608 609 610 611 612 613 614 615 616
        CellClass*const upperCells = new CellClass[nbLevelsAboveRoot+2];

        CellClass*const cellsXAxis = new CellClass[nbLevelsAboveRoot+2];
        CellClass*const cellsYAxis = new CellClass[nbLevelsAboveRoot+2];
        CellClass*const cellsZAxis = new CellClass[nbLevelsAboveRoot+2];
        CellClass*const cellsXYAxis = new CellClass[nbLevelsAboveRoot+2];
        CellClass*const cellsYZAxis = new CellClass[nbLevelsAboveRoot+2];
        CellClass*const cellsXZAxis = new CellClass[nbLevelsAboveRoot+2];

617 618 619 620 621

        // First M2M from level 1 to level 0
        {
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoLeft();
622
            kernels->M2M( &upperCells[nbLevelsAboveRoot+1], octreeIterator.getCurrentBox(), offsetRealTree);
623
        }
624

625 626 627
        // Then M2M from level 0 to level -LIMITE
        {
            CellClass* virtualChild[8];
628 629 630
            for(int idxLevel = nbLevelsAboveRoot ; idxLevel > 0  ; --idxLevel){
                FMemUtils::setall(virtualChild,&upperCells[idxLevel+1],8);
                kernels->M2M( &upperCells[idxLevel], virtualChild, idxLevel + 2);
631
            }
632 633 634 635 636 637 638

            // Cells on the axis of the center should be computed separatly.

            if( usePerDir(DirMinusX) && usePerDir(DirMinusY) ){
                FMemUtils::copyall(cellsZAxis,upperCells,nbLevelsAboveRoot+2);
            }
            else{
639
                 processM2MInIntervals(cellsZAxis,upperCells[nbLevelsAboveRoot+1],
640 641 642 643 644 645
                                    ifDir(DirMinusX,0,1),1,ifDir(DirMinusY,0,1),1,0,1);
            }
            if( usePerDir(DirMinusX) && usePerDir(DirMinusZ) ){
                FMemUtils::copyall(cellsYAxis,upperCells,nbLevelsAboveRoot+2);
            }
            else{
646
                 processM2MInIntervals(cellsYAxis,upperCells[nbLevelsAboveRoot+1],
647 648 649 650 651 652
                                    ifDir(DirMinusX,0,1),1,0,1,ifDir(DirMinusZ,0,1),1);
            }
            if( usePerDir(DirMinusY) && usePerDir(DirMinusZ) ){
                FMemUtils::copyall(cellsXAxis,upperCells,nbLevelsAboveRoot+2);
            }
            else{
653
                 processM2MInIntervals(cellsXAxis,upperCells[nbLevelsAboveRoot+1],
654 655 656 657 658 659
                                    0,1,ifDir(DirMinusY,0,1),1,ifDir(DirMinusZ,0,1),1);
            }

            // Then cells on the spaces should be computed separatly

            if( !usePerDir(DirMinusX) ){
660
                 processM2MInIntervals(cellsYZAxis,upperCells[nbLevelsAboveRoot+1],1,1,0,1,0,1);
661 662 663 664 665
            }
            else {
                FMemUtils::copyall(cellsYZAxis,upperCells,nbLevelsAboveRoot+2);
            }
            if( !usePerDir(DirMinusY) ){
666
                 processM2MInIntervals(cellsXZAxis,upperCells[nbLevelsAboveRoot+1],0,1,1,1,0,1);
667 668 669 670 671
            }
            else {
                FMemUtils::copyall(cellsXZAxis,upperCells,nbLevelsAboveRoot+2);
            }
            if( !usePerDir(DirMinusZ) ){
672
                 processM2MInIntervals(cellsXYAxis,upperCells[nbLevelsAboveRoot+1],0,1,0,1,1,1);
673 674 675 676 677
            }
            else {
                FMemUtils::copyall(cellsXYAxis,upperCells,nbLevelsAboveRoot+2);
            }

678
        }
679

680 681
        // Then M2L at all level
        {
682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714
            CellClass* positionedCells[343];
            memset(positionedCells, 0, 343 * sizeof(CellClass**));

            for(int idxX = ifDir(DirMinusX,-3,0) ; idxX <= ifDir(DirPlusX,3,0) ; ++idxX){
                for(int idxY = ifDir(DirMinusY,-3,0) ; idxY <= ifDir(DirPlusY,3,0) ; ++idxY){
                    for(int idxZ = ifDir(DirMinusZ,-3,0) ; idxZ <= ifDir(DirPlusZ,3,0) ; ++idxZ){
                        if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
                            if(idxX == 0 && idxY == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsZAxis;
                            }
                            else if(idxX == 0 && idxZ == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsYAxis;
                            }
                            else if(idxY == 0 && idxZ == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsXAxis;
                            }
                            else if(idxX == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsYZAxis;
                            }
                            else if(idxY == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsXZAxis;
                            }
                            else if(idxZ == 0){
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = cellsXYAxis;
                            }
                            else{
                                positionedCells[neighIndex(idxX,idxY,idxZ)] = upperCells;
                            }
                        }
                    }
                }
            }

715 716
            // We say that we are in the child index 0
            // So we can compute one time the relative indexes
717 718
            const CellClass* neighbors[343];
            memset(neighbors, 0, sizeof(CellClass*) * 343);
719
            int counter = 0;
720 721 722
            for(int idxX = ifDir(DirMinusX,-3,0) ; idxX <= ifDir(DirPlusX,2,0) ; ++idxX){
                for(int idxY = ifDir(DirMinusY,-3,0) ; idxY <= ifDir(DirPlusY,2,0) ; ++idxY){
                    for(int idxZ = ifDir(DirMinusZ,-3,0) ; idxZ <= ifDir(DirPlusZ,2,0) ; ++idxZ){
723
                        if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
724
                            neighbors[neighIndex(idxX,idxY,idxZ)] = reinterpret_cast<const CellClass*>(~0);
725
                            ++counter;
726 727 728 729
                        }
                    }
                }
            }
730

731
            for(int idxLevel = nbLevelsAboveRoot + 1 ; idxLevel > 1 ; --idxLevel ){
732 733
                for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
                    if(neighbors[idxNeigh]){
734
                        neighbors[idxNeigh] = &positionedCells[idxNeigh][idxLevel];
735
                    }
736
                }
737 738 739 740
                kernels->M2L( &upperCells[idxLevel] , neighbors, counter , idxLevel + 2);
            }

            memset(neighbors, 0, sizeof(CellClass*) * 343);
741 742 743 744
            counter = 0;
            for(int idxX = ifDir(DirMinusX,-2,0) ; idxX <= ifDir(DirPlusX,3,0) ; ++idxX){
                for(int idxY = ifDir(DirMinusY,-2,0) ; idxY <= ifDir(DirPlusY,3,0) ; ++idxY){
                    for(int idxZ = ifDir(DirMinusZ,-2,0) ; idxZ <= ifDir(DirPlusZ,3,0) ; ++idxZ){
745
                        if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
746 747 748
                            const int index = neighIndex(idxX,idxY,idxZ);
                            neighbors[index] = &positionedCells[index][1];
                            ++counter;
749 750 751
                        }
                    }
                }
752
            }
753 754
            kernels->M2L( &upperCells[1] , neighbors, 189, 3);
        }
755

756
        { // compute the border
757 758 759 760 761 762 763
            if( usePerDir(AllDirs) ){
                CellClass leftborder, bottomborder, frontborder, angleborderlb,
                        angleborderfb, angleborderlf, angleborder;

                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                int counter = 0;
764

765 766 767 768 769 770 771 772 773 774 775 776 777 778
                processTopM2MInIntervals( &leftborder, upperCells[nbLevelsAboveRoot+1],     1,1 , 0,1 , 0,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, leftborder,     -2,-2 , -1,1,  -1,1 );
                processTopM2MInIntervals( &bottomborder, upperCells[nbLevelsAboveRoot+1],   0,1 , 0,1 , 1,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, bottomborder,   -1,1  , -1,1,  -2,-2);
                processTopM2MInIntervals( &frontborder, upperCells[nbLevelsAboveRoot+1],    0,1 , 1,1 , 0,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, frontborder,    -1,1  , -2,-2, -1,1 );
                processTopM2MInIntervals( &angleborderlb, upperCells[nbLevelsAboveRoot+1],  1,1 , 0,1 , 1,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, angleborderlb,  -2,-2 , -1,1,  -2,-2);
                processTopM2MInIntervals( &angleborderfb, upperCells[nbLevelsAboveRoot+1],  0,1 , 1,1 , 1,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, angleborderfb,  -1,1 ,  -2,-2, -2,-2);
                processTopM2MInIntervals( &angleborderlf, upperCells[nbLevelsAboveRoot+1],  1,1 , 1,1 , 0,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, angleborderlf,  -2,-2 , -2,-2, -1,1 );
                processTopM2MInIntervals( &angleborder, upperCells[nbLevelsAboveRoot+1],    1,1 , 1,1 , 1,1);
                counter +=  fillM2LVectorFromIntervals(neighbors, angleborder,    -2,-2 , -2,-2, -2,-2);
779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795

                kernels->M2L( &upperCells[0] , neighbors, counter, 2);

                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
                virtualChild[childIndex(0,0,0)] = &upperCells[1];
                kernels->L2L( &upperCells[0], virtualChild, 2);
            }
            else {
                CellClass*const leftborder = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const bottomborder = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const frontborder = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const angleborderlb = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const angleborderfb = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const angleborderlf = new CellClass[nbLevelsAboveRoot+2];
                CellClass*const angleborder = new CellClass[nbLevelsAboveRoot+2];

796 797 798 799 800 801 802
                 processM2MInIntervals( leftborder,   upperCells[nbLevelsAboveRoot+1],     1,1 , 0,1 , 0,1);
                 processM2MInIntervals( bottomborder, upperCells[nbLevelsAboveRoot+1],   0,1 , 0,1 , 1,1);
                 processM2MInIntervals( frontborder,  upperCells[nbLevelsAboveRoot+1],    0,1 , 1,1 , 0,1);
                 processM2MInIntervals( angleborderlb,upperCells[nbLevelsAboveRoot+1],  1,1 , 0,1 , 1,1);
                 processM2MInIntervals( angleborderfb,upperCells[nbLevelsAboveRoot+1],  0,1 , 1,1 , 1,1);
                 processM2MInIntervals( angleborderlf,upperCells[nbLevelsAboveRoot+1],  1,1 , 1,1 , 0,1);
                 processM2MInIntervals( angleborder,  upperCells[nbLevelsAboveRoot+1],    1,1 , 1,1 , 1,1);
803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861

                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                int counter = 0;

                if(usePerDir(DirMinusX) && usePerDir(DirMinusY) && usePerDir(DirMinusZ)){
                    neighbors[neighIndex(-2,-1,-1)] = &leftborder[0];
                    neighbors[neighIndex(-1,-2,-1)] = &frontborder[0];
                    neighbors[neighIndex(-1,-1,-2)] = &bottomborder[0];
                    neighbors[neighIndex(-2,-2,-1)] = &angleborderlf[0];
                    neighbors[neighIndex(-2,-1,-2)] = &angleborderlb[0];
                    neighbors[neighIndex(-1,-2,-2)] = &angleborderfb[0];
                    neighbors[neighIndex(-2,-2,-2)] = &angleborder[0];
                    counter += 7;
                }
                if(usePerDir(DirMinusX) && usePerDir(DirPlusY) && usePerDir(DirMinusZ)){
                    neighbors[neighIndex(-2,1,-1)] = &leftborder[0];
                    neighbors[neighIndex(-1,1,-2)] = &bottomborder[0];
                    neighbors[neighIndex(-2,1,-2)] = &angleborderlb[0];
                    counter += 3;
                }
                if(usePerDir(DirMinusX) && usePerDir(DirMinusY) && usePerDir(DirPlusZ)){
                    neighbors[neighIndex(-2,-1,1)] = &leftborder[0];
                    neighbors[neighIndex(-1,-2,1)] = &frontborder[0];
                    neighbors[neighIndex(-2,-2,1)] = &angleborderlf[0];
                    counter += 3;
                }
                if(usePerDir(DirMinusX) && usePerDir(DirPlusY) && usePerDir(DirPlusZ)){
                    neighbors[neighIndex(-2,1,1)] = &leftborder[0];
                    counter += 1;
                }
                if(usePerDir(DirPlusX) && usePerDir(DirMinusY) && usePerDir(DirMinusZ)){
                    neighbors[neighIndex(1,-2,-1)] = &frontborder[0];
                    neighbors[neighIndex(1,-1,-2)] = &bottomborder[0];
                    neighbors[neighIndex(1,-2,-2)] = &angleborderfb[0];
                    counter += 3;
                }
                if(usePerDir(DirPlusX) && usePerDir(DirMinusY) && usePerDir(DirPlusZ)){
                    neighbors[neighIndex(1,-2,1)] = &frontborder[0];
                    counter += 1;
                }
                if(usePerDir(DirPlusX) && usePerDir(DirPlusY) && usePerDir(DirMinusZ)){
                    neighbors[neighIndex(1,1,-2)] = &bottomborder[0];
                    counter += 1;
                }

                CellClass centerXFace;
                CellClass centerYFace;
                CellClass centerZFace;
                CellClass centerXZFace;
                CellClass centerXYFace;
                CellClass centerYXFace;
                CellClass centerYZFace;
                CellClass centerZXFace;
                CellClass centerZYFace;
                CellClass angleXZFace;
                CellClass angleXYFace;
                CellClass angleYZFace;

862
                if(usePerDir(DirMinusX)){
863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925
                    {
                        CellClass* virtualChild[8];
                        memset(virtualChild, 0, 8*sizeof(CellClass*));
                        if(usePerDir(DirPlusY) && usePerDir(DirPlusZ))  virtualChild[childIndex(1,1,1)] = &leftborder[1];
                        if(usePerDir(DirMinusY) && usePerDir(DirPlusZ)) virtualChild[childIndex(1,0,1)] = &leftborder[1];
                        else if( usePerDir(DirPlusZ))                   virtualChild[childIndex(1,0,1)] = &angleborderlf[1];
                        if(usePerDir(DirPlusY) && usePerDir(DirMinusZ)) virtualChild[childIndex(1,1,0)] = &leftborder[1];
                        else if(usePerDir(DirPlusY) )                   virtualChild[childIndex(1,1,0)] = &angleborderlb[1];

                        if(usePerDir(DirMinusY) && usePerDir(DirMinusZ)) virtualChild[childIndex(1,0,0)] = &leftborder[1];
                        else if(usePerDir(DirMinusZ))                    virtualChild[childIndex(1,0,0)] = &angleborderlf[1];
                        else if(usePerDir(DirMinusY))                    virtualChild[childIndex(1,0,0)] = &angleborderlb[1];
                        else                                             virtualChild[childIndex(1,0,0)] = &angleborder[1];

                        kernels->M2M( &centerXFace, virtualChild, 2);
                        neighbors[neighIndex(-2,0,0)] = &centerXFace;
                        counter += 1;
                    }
                    if(usePerDir(DirMinusZ) || usePerDir(DirPlusZ)){
                        if(usePerDir(DirY)){
                            centerXZFace = leftborder[0];
                        }
                        else{
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusY)){
                                virtualChild[childIndex(1,1,0)] = &leftborder[1];
                                virtualChild[childIndex(1,1,1)] = &leftborder[1];
                            }
                            if(usePerDir(DirMinusY)){
                                virtualChild[childIndex(1,0,0)] = &leftborder[1];
                                virtualChild[childIndex(1,0,1)] = &leftborder[1];
                            }
                            else{
                                virtualChild[childIndex(1,0,0)] = &angleborderlf[1];
                                virtualChild[childIndex(1,0,1)] = &angleborderlf[1];
                            }
                            kernels->M2M( &centerXZFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusZ) ){
                            neighbors[neighIndex(-2,0,1)] = &centerXZFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusZ) ){
                            neighbors[neighIndex(-2,0,-1)] = &centerXZFace;
                            counter += 1;

                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusY)){
                                virtualChild[childIndex(1,1,1)] = &angleborderlb[1];
                            }
                            if(usePerDir(DirMinusY)){
                                virtualChild[childIndex(1,0,1)] = &angleborderlb[1];
                            }
                            else{
                                virtualChild[childIndex(1,0,1)] = &angleborder[1];
                            }
                            kernels->M2M( &angleXZFace, virtualChild, 2);

                            neighbors[neighIndex(-2,0,-2)] = &angleXZFace;
                            counter += 1;
                        }
926
                    }
927 928 929
                    if(usePerDir(DirMinusY) || usePerDir(DirPlusY)){
                        if(usePerDir(DirZ)){
                            centerXYFace = leftborder[0];
930 931
                        }
                        else{
932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusZ)){
                                virtualChild[childIndex(1,0,1)] = &leftborder[1];
                                virtualChild[childIndex(1,1,1)] = &leftborder[1];
                            }
                            if(usePerDir(DirMinusZ)){
                                virtualChild[childIndex(1,0,0)] = &leftborder[1];
                                virtualChild[childIndex(1,1,0)] = &leftborder[1];
                            }
                            else{
                                virtualChild[childIndex(1,0,0)] = &angleborderlb[1];
                                virtualChild[childIndex(1,1,0)] = &angleborderlb[1];
                            }
                            kernels->M2M( &centerXYFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusY) ){
                            neighbors[neighIndex(-2,1,0)] = &centerXYFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusY) ){
                            neighbors[neighIndex(-2,-1,0)] = &centerXYFace;
                            counter += 1;

                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusZ)){
                                virtualChild[childIndex(1,1,1)] = &angleborderlf[1];
                            }
                            if(usePerDir(DirMinusZ)){
                                virtualChild[childIndex(1,1,0)] = &angleborderlf[1];
                            }
                            else{
                                virtualChild[childIndex(1,1,0)] = &angleborder[1];
                            }
                            kernels->M2M( &angleXYFace, virtualChild, 2);

                            neighbors[neighIndex(-2,-2,0)] = &angleXYFace;
                            counter += 1;
971 972 973
                        }
                    }
                }
974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991
                if(usePerDir(DirMinusY)){
                    {
                        CellClass* virtualChild[8];
                        memset(virtualChild, 0, 8*sizeof(CellClass*));
                        if(usePerDir(DirPlusX) && usePerDir(DirPlusZ))  virtualChild[childIndex(1,1,1)] = &frontborder[1];
                        if(usePerDir(DirMinusX) && usePerDir(DirPlusZ)) virtualChild[childIndex(0,1,1)] = &frontborder[1];
                        else if(usePerDir(DirPlusZ))                    virtualChild[childIndex(0,1,1)] = &angleborderlf[1];
                        if(usePerDir(DirPlusX) && usePerDir(DirMinusZ)) virtualChild[childIndex(1,1,0)] = &frontborder[1];
                        else if(usePerDir(DirPlusX))                    virtualChild[childIndex(1,1,0)] = &angleborderfb[1];

                        if(usePerDir(DirMinusX) && usePerDir(DirMinusZ)) virtualChild[childIndex(0,1,0)] = &frontborder[1];
                        else if(usePerDir(DirMinusZ))                    virtualChild[childIndex(0,1,0)] = &angleborderlf[1];
                        else if(usePerDir(DirMinusX))                    virtualChild[childIndex(0,1,0)] = &angleborderfb[1];
                        else                                             virtualChild[childIndex(0,1,0)] = &angleborder[1];

                        kernels->M2M( &centerYFace, virtualChild, 2);
                        neighbors[neighIndex(0,-2,0)] = &centerYFace;
                        counter += 1;
992
                    }
993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067
                    if(usePerDir(DirMinusZ) || usePerDir(DirPlusZ)){
                        if(usePerDir(DirX)){
                            centerYZFace = frontborder[0];
                        }
                        else{
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusX)){
                                virtualChild[childIndex(1,1,0)] = &frontborder[1];
                                virtualChild[childIndex(1,1,1)] = &frontborder[1];
                            }
                            if(usePerDir(DirMinusX)){
                                virtualChild[childIndex(0,1,0)] = &frontborder[1];
                                virtualChild[childIndex(0,1,1)] = &frontborder[1];
                            }
                            else{
                                virtualChild[childIndex(0,1,0)] = &angleborderlf[1];
                                virtualChild[childIndex(0,1,1)] = &angleborderlf[1];
                            }
                            kernels->M2M( &centerYZFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusZ) ){
                            neighbors[neighIndex(0,-2,1)] = &centerYZFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusZ) ){
                            neighbors[neighIndex(0,-2,-1)] = &centerYZFace;
                            counter += 1;

                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusX)){
                                virtualChild[childIndex(1,1,1)] = &angleborderfb[1];
                            }
                            if(usePerDir(DirMinusX)){
                                virtualChild[childIndex(0,1,1)] = &angleborderfb[1];
                            }
                            else{
                                virtualChild[childIndex(0,1,1)] = &angleborder[1];
                            }
                            kernels->M2M( &angleYZFace, virtualChild, 2);

                            neighbors[neighIndex(0,-2,-2)] = &angleYZFace;
                            counter += 1;
                        }
                    }
                    if(usePerDir(DirMinusX) || usePerDir(DirPlusX)){
                        if(usePerDir(DirZ)){
                            centerYXFace = frontborder[0];
                        }
                        else{
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusZ)){
                                virtualChild[childIndex(0,1,1)] = &frontborder[1];
                                virtualChild[childIndex(1,1,1)] = &frontborder[1];
                            }
                            if(usePerDir(DirMinusZ)){
                                virtualChild[childIndex(0,1,0)] = &frontborder[1];
                                virtualChild[childIndex(1,1,0)] = &frontborder[1];
                            }
                            else{
                                virtualChild[childIndex(0,1,0)] = &angleborderfb[1];
                                virtualChild[childIndex(1,1,0)] = &angleborderfb[1];
                            }
                            kernels->M2M( &centerYXFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusX) ){
                            neighbors[neighIndex(1,-2,0)] = &centerYXFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusX) ){
                            neighbors[neighIndex(-1,-2,0)] = &centerYXFace;
                            counter += 1;
                        }
1068 1069
                    }
                }
1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148
                if(usePerDir(DirMinusZ)){
                    {
                        CellClass* virtualChild[8];
                        memset(virtualChild, 0, 8*sizeof(CellClass*));
                        if(usePerDir(DirPlusX) && usePerDir(DirPlusY))  virtualChild[childIndex(1,1,1)] = &bottomborder[1];
                        if(usePerDir(DirMinusX) && usePerDir(DirPlusY)) virtualChild[childIndex(0,1,1)] = &bottomborder[1];
                        else if( usePerDir(DirPlusY))                   virtualChild[childIndex(0,1,1)] = &angleborderlb[1];
                        if(usePerDir(DirPlusX) && usePerDir(DirMinusY)) virtualChild[childIndex(1,0,1)] = &bottomborder[1];
                        else if(usePerDir(DirPlusX) )                   virtualChild[childIndex(1,0,1)] = &angleborderfb[1];

                        if(usePerDir(DirMinusX) && usePerDir(DirMinusY)) virtualChild[childIndex(0,0,1)] = &bottomborder[1];
                        else if(usePerDir(DirMinusY))                    virtualChild[childIndex(0,0,1)] = &angleborderfb[1];
                        else if(usePerDir(DirMinusX))                    virtualChild[childIndex(0,0,1)] = &angleborderlb[1];
                        else                                             virtualChild[childIndex(0,0,1)] = &angleborder[1];

                        kernels->M2M( &centerZFace, virtualChild, 2);
                        neighbors[neighIndex(0,0,-2)] = &centerZFace;
                        counter += 1;
                    }
                    if(usePerDir(DirMinusY) || usePerDir(DirPlusY)){
                        if(usePerDir(DirX)){
                            centerZYFace = bottomborder[0];
                        }
                        else{
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusX)){
                                virtualChild[childIndex(1,0,1)] = &bottomborder[1];
                                virtualChild[childIndex(1,1,1)] = &bottomborder[1];
                            }
                            if(usePerDir(DirMinusX)){
                                virtualChild[childIndex(0,0,1)] = &bottomborder[1];
                                virtualChild[childIndex(0,1,1)] = &bottomborder[1];
                            }
                            else{
                                virtualChild[childIndex(0,0,1)] = &angleborderlb[1];
                                virtualChild[childIndex(0,1,1)] = &angleborderlb[1];
                            }
                            kernels->M2M( &centerZYFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusY) ){
                            neighbors[neighIndex(0,1,-2)] = &centerZYFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusY) ){
                            neighbors[neighIndex(0,-1,-2)] = &centerZYFace;
                            counter += 1;
                        }
                    }
                    if(usePerDir(DirMinusX) || usePerDir(DirPlusX)){
                        if(usePerDir(DirY)){
                            centerZXFace = bottomborder[0];
                        }
                        else{
                            CellClass* virtualChild[8];
                            memset(virtualChild, 0, 8*sizeof(CellClass*));
                            if(usePerDir(DirPlusY)){
                                virtualChild[childIndex(0,1,1)] = &bottomborder[1];
                                virtualChild[childIndex(1,1,1)] = &bottomborder[1];
                            }
                            if(usePerDir(DirMinusY)){
                                virtualChild[childIndex(0,0,1)] = &bottomborder[1];
                                virtualChild[childIndex(1,0,1)] = &bottomborder[1];
                            }
                            else{
                                virtualChild[childIndex(0,0,1)] = &angleborderlf[1];
                                virtualChild[childIndex(0,1,1)] = &angleborderlf[1];
                            }
                            kernels->M2M( &centerZXFace, virtualChild, 2);
                        }
                        if( usePerDir(DirPlusX) ){
                            neighbors[neighIndex(1,0,-2)] = &centerZXFace;
                            counter += 1;
                        }
                        if( usePerDir(DirMinusX) ){
                            neighbors[neighIndex(-1,0,-2)] = &centerZXFace;
                            counter += 1;
                        }
                    }
1149 1150
                }

1151
                // M2L for border
1152 1153
                kernels->M2L( &upperCells[0] , neighbors, counter, 2);

1154
                // L2L from border M2L to top of tree
1155 1156 1157 1158 1159
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
                virtualChild[childIndex(0,0,0)] = &upperCells[1];
                kernels->L2L( &upperCells[0], virtualChild, 2);

1160
                // dealloc
1161 1162 1163 1164 1165 1166 1167 1168
                delete[] leftborder;
                delete[] bottomborder;
                delete[] frontborder;
                delete[] angleborderlb;
                delete[] angleborderfb;
                delete[] angleborderlf;
                delete[] angleborder;
            }
1169

1170 1171 1172 1173 1174 1175
        }

        // Finally L2L until level 0
        {
            CellClass* virtualChild[8];
            memset(virtualChild, 0, sizeof(CellClass*) * 8);
1176
            for(int idxLevel = 1 ; idxLevel <= nbLevelsAboveRoot  ; ++idxLevel){
1177
                virtualChild[childIndex(1,1,1)] = &upperCells[idxLevel+1];
1178
                kernels->L2L( &upperCells[idxLevel], virtualChild, idxLevel + 2);
1179 1180 1181 1182 1183 1184 1185
            }
        }

        // L2L from 0 to level 1
        {
            typename OctreeClass::Iterator octreeIterator(tree);
            octreeIterator.gotoLeft();
1186
            kernels->L2L( &upperCells[nbLevelsAboveRoot+1], octreeIterator.getCurrentBox(), offsetRealTree);
1187
        }
1188

1189
        delete[] upperCells;
1190

1191 1192 1193 1194 1195 1196
        delete[] cellsXAxis;
        delete[] cellsYAxis;
        delete[] cellsZAxis;
        delete[] cellsXYAxis;
        delete[] cellsYZAxis;
        delete[] cellsXZAxis;
1197 1198

        FDEBUG( FDebug::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
1199 1200 1201
    }


1202 1203 1204 1205
};


#endif // FFMMALGORITHMPERIODIC_HPP