FFmmAlgorithmPeriodic.hpp 24.5 KB
Newer Older
1
// See LICENCE file at project root
2 3
#ifndef FFMMALGORITHMPERIODIC_HPP
#define FFMMALGORITHMPERIODIC_HPP
4

5 6

#include "../Utils/FGlobal.hpp"
7
#include "../Utils/FGlobalPeriodic.hpp"
8
#include "../Utils/FAssert.hpp"
9
#include "../Utils/FLog.hpp"
10

11
#include "../Utils/FTic.hpp"
12
#include "../Utils/FMemUtils.hpp"
13 14 15 16

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

17
#include "FCoreCommon.hpp"
18 19 20 21 22 23 24

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFmmAlgorithmPeriodic
* @brief
* Please read the license
*
25
* This class is a basic FMM algorithm with periodic behavior
26 27
* It just iterates on a tree and call the kernels with good arguments.
*
COULAUD Olivier's avatar
COULAUD Olivier committed
28
* Of course this class does not deallocate pointer given in arguments.
29
*/
30
template<class FReal, class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
31
class FFmmAlgorithmPeriodic : public FAbstractAlgorithm{
32

33 34
    OctreeClass* const tree;        //< The octree to work on
    KernelClass* kernels;           //< The kernels
35

COULAUD Olivier's avatar
COULAUD Olivier committed
36
    const int OctreeHeight;         //< The height of the octree (real height)
37 38 39
    const int nbLevelsAboveRoot;    //< The nb of level the user ask to go above the tree (>= -1)
    const int offsetRealTree;       //< nbLevelsAboveRoot GetFackLevel

40
    const int leafLevelSeperationCriteria;
41 42 43 44 45 46

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
COULAUD Olivier's avatar
COULAUD Olivier committed
47
      * @param inUpperLevel this parameter defines the behavior of the periodicity refer to the main doc
48
      *
49
      */
50
    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inLeafLevelSeperationCriteria = 1)
51
        : tree(inTree) , kernels(nullptr), OctreeHeight(tree->getHeight()),
52
          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 2), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
53

54 55
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
56
        FAssertLF(leafLevelSeperationCriteria < 3, "Separation criteria should be < 3");
57

58 59
        FAbstractAlgorithm::setNbLevelsInTree(extendedTreeHeight());

60
        FLOG(FLog::Controller << "FFmmAlgorithmPeriodic\n");
61 62 63 64 65 66
    }

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

67 68 69 70
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

71 72 73 74 75 76 77 78



    long long int theoricalRepetition() const {
        if( nbLevelsAboveRoot == -1 ){
            // we know it is 3 (-1;+1)
            return 3;
        }
79
        return 6 * (1 << nbLevelsAboveRoot);
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
    }


    void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
        if( nbLevelsAboveRoot == -1 ){
            // We know it is (-1;1)
            min->setPosition(-1,-1,-1);
            max->setPosition(1,1,1);
        }
        else{
            const int halfRepeated = int(theoricalRepetition()/2);
            min->setPosition(-halfRepeated,-halfRepeated,-halfRepeated);
            // if we repeat the box 8 times, we go from [-4 to 3]
            max->setPosition(halfRepeated-1,halfRepeated-1,halfRepeated-1);
        }
    }


    FReal extendedBoxWidth() const {
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
        if( nbLevelsAboveRoot == -1 ){
            return tree->getBoxWidth()*2;
        }
        else{
            return tree->getBoxWidth() * FReal(4<<(nbLevelsAboveRoot));
        }
    }

    FReal extendedBoxWidthBoundary() const {
        if( nbLevelsAboveRoot == -1 ){
            return tree->getBoxWidth()*4;
        }
        else{
            return tree->getBoxWidth() * FReal(8<<(nbLevelsAboveRoot));
        }
114 115 116 117 118 119 120 121 122
    }

    /** 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
      */
123
    FPoint<FReal> extendedBoxCenter() const {
124 125 126 127 128 129 130 131 132 133 134 135
        if( nbLevelsAboveRoot == -1 ){
            const FReal originalBoxWidth            = tree->getBoxWidth();
            const FPoint<FReal> originalBoxCenter   = tree->getBoxCenter();
            const FReal originalBoxWidthDiv2        = originalBoxWidth/2.0;
            return FPoint<FReal>( originalBoxCenter.getX() + originalBoxWidthDiv2,
                                         originalBoxCenter.getY() + originalBoxWidthDiv2,
                                         originalBoxCenter.getZ() + originalBoxWidthDiv2);
        }
        else{
            const FReal originalBoxWidth     = tree->getBoxWidth();
            const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
            const FPoint<FReal> originalBoxCenter   = tree->getBoxCenter();
136

137 138
            const FReal offset = extendedBoxWidth()/FReal(2.0);
            return FPoint<FReal>( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
139 140
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
        }
    }

    FPoint<FReal> extendedBoxCenterBoundary() const {
        if( nbLevelsAboveRoot == -1 ){
            const FReal originalBoxWidth            = tree->getBoxWidth();
            const FPoint<FReal> originalBoxCenter   = tree->getBoxCenter();
            const FReal originalBoxWidthDiv2        = originalBoxWidth/2.0;
            return FPoint<FReal>( originalBoxCenter.getX() + originalBoxWidthDiv2,
                                         originalBoxCenter.getY() + originalBoxWidthDiv2,
                                         originalBoxCenter.getZ() + originalBoxWidthDiv2);
        }
        else{
            const FReal originalBoxWidth     = tree->getBoxWidth();
            const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
            const FPoint<FReal> originalBoxCenter   = tree->getBoxCenter();

            return FPoint<FReal>( originalBoxCenter.getX() + originalBoxWidthDiv2,
                       originalBoxCenter.getY() + originalBoxWidthDiv2,
                       originalBoxCenter.getZ() + originalBoxWidthDiv2);
        }
162 163 164 165 166 167 168 169 170 171 172 173 174
    }

    /** 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
      */
    int extendedTreeHeight() const {
        // The real height
        return OctreeHeight + offsetRealTree;
    }

175 176 177 178 179
    int extendedTreeHeightBoundary() const {
        // The real height
        return OctreeHeight + offsetRealTree + 1;
    }

180
protected:
181 182 183 184
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
185
    void executeCore(const unsigned operationsToProceed) override {
186
        FAssertLF(kernels, "kernels cannot be null");
187

188
        if(operationsToProceed & FFmmP2M) bottomPass();
189

190
        if(operationsToProceed & FFmmM2M) upwardPass();
191

192 193 194 195 196
        if(operationsToProceed & FFmmM2L){
            transferPass();
            // before downward pass we have to perform the periodicity
            processPeriodicLevels();
        }
197

198
        if(operationsToProceed & FFmmL2L) downardPass();
199

200
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass((operationsToProceed & FFmmP2P),(operationsToProceed & FFmmL2P));
201 202
    }

203

204 205 206 207 208 209
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
210
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
211 212
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
213 214 215 216 217 218 219 220

        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
221
            FLOG(computationCounter.tic());
222
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
223
            FLOG(computationCounter.tac());
224 225
        } while(octreeIterator.moveRight());

226
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << " s)\n" );
227
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
228 229 230 231 232 233 234 235
    }

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

    /** M2M */
    void upwardPass(){
236
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
237 238
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
239 240 241 242 243 244 245 246 247 248

        // 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 ){
249
            FLOG(FTic counterTimeLevel);
250
            const int fackLevel = idxLevel + offsetRealTree;
251 252 253 254
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
255
                FLOG(computationCounter.tic());
256
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
257
                FLOG(computationCounter.tac());
258 259 260 261
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
262
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
263 264 265
        }


266
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << " s)\n" );
267
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
268 269 270
    }

    /////////////////////////////////////////////////////////////////////////////
271
    // Transfer
272 273 274
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
275
    void transferPass(){
276
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
277 278
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
279 280 281 282

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

283 284
        const CellClass* neighbors[342];
        int neighborPositions[342];
285 286 287

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
288
            FLOG(FTic counterTimeLevel);
289
            const int fackLevel = idxLevel + offsetRealTree;
290
            const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
291 292
            // for each cells
            do{
293
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
294
                FLOG(computationCounter.tic());
295
                if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, neighborPositions, counter, fackLevel);
296
                FLOG(computationCounter.tac());
297 298 299
            } while(octreeIterator.moveRight());
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
300

301
            FLOG(computationCounter.tic());
302
            kernels->finishedLevelM2L(fackLevel);
303
            FLOG(computationCounter.tac());
304
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
305
        }
306
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << " s)\n" );
307
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
308
    }
309

310 311 312 313 314 315
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
316
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
317 318
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
319 320 321

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

323 324 325
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
326
            FLOG(FTic counterTimeLevel);
327 328
            const int fackLevel = idxLevel + offsetRealTree;

329 330
            // for each cells
            do{
331
                FLOG(computationCounter.tic());
332
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
333
                FLOG(computationCounter.tac());
334 335 336 337
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
338
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
339 340
        }

341
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (L2L) = "  << counterTime.tacAndElapsed() << " s)\n" );
342
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
343

344 345 346 347 348 349 350 351

    }

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

    /** P2P */
352
    void directPass(const bool p2pEnabled, const bool l2pEnabled){
353
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
354 355 356
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
357 358

        const int heightMinusOne = OctreeHeight - 1;
359
        const FReal boxWidth = tree->getBoxWidth();
360 361 362 363

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
364 365 366
        ContainerClass* neighbors[26];
        FTreeCoordinate offsets[26];
        int neighborPositions[26];
367
        bool hasPeriodicLeaves;
368 369
        // for each leafs
        do{
370 371 372 373 374 375 376 377
            if(l2pEnabled){
                FLOG(computationCounterL2P.tic());
                kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
                FLOG(computationCounterL2P.tac());
            }
            if(p2pEnabled){
                // need the current particles and neighbors particles
                const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
378
                const int counter = tree->getPeriodicLeafsNeighbors( neighbors, neighborPositions, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
379 380 381
                int periodicNeighborsCounter = 0;

                if(hasPeriodicLeaves){
382 383
                    ContainerClass* periodicNeighbors[26];
                    int periodicNeighborPositions[26];
384

385 386
                    for(int idxNeig = 0 ; idxNeig < counter ; ++idxNeig){
                        if( !offsets[idxNeig].equals(0,0,0) ){
387
                            FAssertLF(octreeIterator.getCurrentListTargets() != neighbors[idxNeig]);
388
                            // Put periodic neighbors into other array
BRAMAS Berenger's avatar
BRAMAS Berenger committed
389 390 391
                            FReal*const positionsX = neighbors[idxNeig]->getWPositions()[0];
                            FReal*const positionsY = neighbors[idxNeig]->getWPositions()[1];
                            FReal*const positionsZ = neighbors[idxNeig]->getWPositions()[2];
392

BRAMAS Berenger's avatar
BRAMAS Berenger committed
393
                            for(FSize idxPart = 0; idxPart < neighbors[idxNeig]->getNbParticles() ; ++idxPart){
394 395 396 397
                                positionsX[idxPart] += boxWidth * FReal(offsets[idxNeig].getX());
                                positionsY[idxPart] += boxWidth * FReal(offsets[idxNeig].getY());
                                positionsZ[idxPart] += boxWidth * FReal(offsets[idxNeig].getZ());
                            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
398

399
                            offsets[periodicNeighborsCounter] = offsets[idxNeig];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
400 401 402
                            periodicNeighbors[periodicNeighborsCounter] = neighbors[idxNeig];
                            periodicNeighborPositions[periodicNeighborsCounter] = neighborPositions[idxNeig];
                            ++periodicNeighborsCounter;
403
                        }
404 405 406 407
                        else{
                            neighbors[idxNeig-periodicNeighborsCounter] = neighbors[idxNeig];
                            neighborPositions[idxNeig-periodicNeighborsCounter] = neighborPositions[idxNeig];
                        }
408
                    }
409

410 411
                    FLOG(computationCounterP2P.tic());
                    kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
412
                                 octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborPositions, periodicNeighborsCounter);
413 414
                    FLOG(computationCounterP2P.tac());

415 416 417 418
                    for(int idxNeig = 0 ; idxNeig < periodicNeighborsCounter ; ++idxNeig){
                        FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
                        FReal*const positionsY = periodicNeighbors[idxNeig]->getWPositions()[1];
                        FReal*const positionsZ = periodicNeighbors[idxNeig]->getWPositions()[2];
419

420 421 422 423
                        for(FSize 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());
424 425 426
                        }
                    }
                }
427

428 429
                FLOG(computationCounterP2P.tic());
                kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
430
                             octreeIterator.getCurrentListSrc(), neighbors, neighborPositions, counter - periodicNeighborsCounter);
431 432
                FLOG(computationCounterP2P.tac());
            }
433 434 435
        } while(octreeIterator.moveRight());


436
        FLOG( FLog::Controller << "\tFinished (@Direct Pass (L2P + P2P) = "  << counterTime.tacAndElapsed() << " s)\n" );
437 438
        FLOG( FLog::Controller << "\t\t Computation L2P : " << computationCounterL2P.cumulated() << " s\n" );
        FLOG( FLog::Controller << "\t\t Computation P2P : " << computationCounterP2P.cumulated() << " s\n" );
439 440 441 442 443 444 445

    }

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

446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
    /** 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);
    }

    /** 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
      */
467
    void processPeriodicLevels(){
468
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
469
        FLOG(FTic counterTime);
470

BRAMAS Berenger's avatar
BRAMAS Berenger committed
471 472 473 474 475 476 477 478
        if( nbLevelsAboveRoot != -1 ){
            // we will use offsetRealTree-1 cells but for simplicity allocate offsetRealTree
            // upperCells[offsetRealTree-1] is root cell
            CellClass*const upperCells = new CellClass[offsetRealTree];
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->M2M( &upperCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
479
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
480 481 482 483 484
            {
                CellClass* virtualChild[8];
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
                    FMemUtils::setall(virtualChild,&upperCells[idxLevel],8);
                    kernels->M2M( &upperCells[idxLevel-1], virtualChild, idxLevel);
485
                }
486
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
487
            CellClass*const downerCells = new CellClass[offsetRealTree];
488

BRAMAS Berenger's avatar
BRAMAS Berenger committed
489 490
            {
                const int idxUpperLevel = 2;
491

492 493
                const CellClass* neighbors[342];
                int neighborPositions[342];
494
                int counter = 0;
495 496 497
                for(int idxX = -3 ; idxX <= 2 ; ++idxX){
                    for(int idxY = -3 ; idxY <= 2 ; ++idxY){
                        for(int idxZ = -3 ; idxZ <= 2 ; ++idxZ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
498
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
499 500
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
501
                                ++counter;
502
                            }
503 504 505
                        }
                    }
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
506
                // compute M2L
507
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, neighborPositions, counter, idxUpperLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
508
            }
509

BRAMAS Berenger's avatar
BRAMAS Berenger committed
510
            for(int idxUpperLevel = 3 ; idxUpperLevel <= offsetRealTree ; ++idxUpperLevel){
511 512
                const CellClass* neighbors[342];
                int neighborPositions[342];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
513 514 515 516 517
                int counter = 0;
                for(int idxX = -2 ; idxX <= 3 ; ++idxX){
                    for(int idxY = -2 ; idxY <= 3 ; ++idxY){
                        for(int idxZ = -2 ; idxZ <= 3 ; ++idxZ){
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
518 519
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
520
                                ++counter;
521 522 523
                            }
                        }
                    }
524 525
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
526
                // compute M2L
527
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, neighborPositions, counter, idxUpperLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
528
            }
529

BRAMAS Berenger's avatar
BRAMAS Berenger committed
530
            {
531 532
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
533
                for(int idxLevel = 2 ; idxLevel < offsetRealTree-1  ; ++idxLevel){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
534 535 536
                    virtualChild[0] = &downerCells[idxLevel];
                    kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
                }
537
            }
538

539 540 541 542 543 544 545 546
            {
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
                const int idxLevel = offsetRealTree-1;
                virtualChild[7] = &downerCells[idxLevel];
                kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
547 548 549 550 551
            // L2L from 0 to level 1
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->L2L( &downerCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
552 553
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
554 555
            delete[] upperCells;
            delete[] downerCells;
556
        }
557

558
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << " s)\n" );
559 560 561
    }


562 563 564 565
};


#endif // FFMMALGORITHMPERIODIC_HPP