FFmmAlgorithmPeriodic.hpp 25.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
// ===================================================================================
// 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".
// ===================================================================================
16 17
#ifndef FFMMALGORITHMPERIODIC_HPP
#define FFMMALGORITHMPERIODIC_HPP
18

19 20

#include "../Utils/FGlobal.hpp"
21
#include "../Utils/FGlobalPeriodic.hpp"
22
#include "../Utils/FAssert.hpp"
23
#include "../Utils/FLog.hpp"
24

25
#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
* It just iterates on a tree and call the kernels with good arguments.
*
COULAUD Olivier's avatar
COULAUD Olivier committed
42
* Of course this class does not deallocate pointer given in arguments.
43
*/
44
template<class FReal, class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
45
class FFmmAlgorithmPeriodic : public FAbstractAlgorithm{
46

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

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

54
    const int leafLevelSeperationCriteria;
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
COULAUD Olivier's avatar
COULAUD Olivier committed
61
      * @param inUpperLevel this parameter defines the behavior of the periodicity refer to the main doc
62
      *
63
      */
64
    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inLeafLevelSeperationCriteria = 1)
65
        : tree(inTree) , kernels(nullptr), OctreeHeight(tree->getHeight()),
66
          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 2), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
67

68 69
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
70
        FAssertLF(leafLevelSeperationCriteria < 3, "Separation criteria should be < 3");
71

72 73
        FAbstractAlgorithm::setNbLevelsInTree(extendedTreeHeight());

74
        FLOG(FLog::Controller << "FFmmAlgorithmPeriodic\n");
75 76 77 78 79 80
    }

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

81 82 83 84
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

85 86 87 88 89 90 91 92



    long long int theoricalRepetition() const {
        if( nbLevelsAboveRoot == -1 ){
            // we know it is 3 (-1;+1)
            return 3;
        }
93
        return 6 * (1 << nbLevelsAboveRoot);
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
    }


    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 {
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
        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));
        }
128 129 130 131 132 133 134 135 136
    }

    /** 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
      */
137
    FPoint<FReal> extendedBoxCenter() const {
138 139 140 141 142 143 144 145 146 147 148 149
        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();
150

151 152
            const FReal offset = extendedBoxWidth()/FReal(2.0);
            return FPoint<FReal>( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
153 154
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
        }
    }

    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);
        }
176 177 178 179 180 181 182 183 184 185 186 187 188
    }

    /** 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;
    }

189 190 191 192 193
    int extendedTreeHeightBoundary() const {
        // The real height
        return OctreeHeight + offsetRealTree + 1;
    }

194
protected:
195 196 197 198
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
199
    void executeCore(const unsigned operationsToProceed) override {
200
        FAssertLF(kernels, "kernels cannot be null");
201

202
        if(operationsToProceed & FFmmP2M) bottomPass();
203

204
        if(operationsToProceed & FFmmM2M) upwardPass();
205

206 207 208 209 210
        if(operationsToProceed & FFmmM2L){
            transferPass();
            // before downward pass we have to perform the periodicity
            processPeriodicLevels();
        }
211

212
        if(operationsToProceed & FFmmL2L) downardPass();
213

214
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass((operationsToProceed & FFmmP2P),(operationsToProceed & FFmmL2P));
215 216
    }

217

218 219 220 221 222 223
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
224
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
225 226
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
227 228 229 230 231 232 233 234

        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
235
            FLOG(computationCounter.tic());
236
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
237
            FLOG(computationCounter.tac());
238 239
        } while(octreeIterator.moveRight());

240 241
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
242 243 244 245 246 247 248 249
    }

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

    /** M2M */
    void upwardPass(){
250
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
251 252
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
253 254 255 256 257 258 259 260 261 262

        // 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 ){
263
            FLOG(FTic counterTimeLevel);
264
            const int fackLevel = idxLevel + offsetRealTree;
265 266 267 268
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
269
                FLOG(computationCounter.tic());
270
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
271
                FLOG(computationCounter.tac());
272 273 274 275
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
276
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
277 278 279
        }


280 281
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
282 283 284
    }

    /////////////////////////////////////////////////////////////////////////////
285
    // Transfer
286 287 288
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
289
    void transferPass(){
290
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
291 292
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
293 294 295 296

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

297 298
        const CellClass* neighbors[342];
        int neighborPositions[342];
299 300 301

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
302
            FLOG(FTic counterTimeLevel);
303
            const int fackLevel = idxLevel + offsetRealTree;
304
            const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
305 306
            // for each cells
            do{
307
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
308
                FLOG(computationCounter.tic());
309
                if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, neighborPositions, counter, fackLevel);
310
                FLOG(computationCounter.tac());
311 312 313
            } while(octreeIterator.moveRight());
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
314

315
            FLOG(computationCounter.tic());
316
            kernels->finishedLevelM2L(fackLevel);
317
            FLOG(computationCounter.tac());
318
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
319
        }
320 321
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
322
    }
323

324 325 326 327 328 329
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
330
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
331 332
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
333 334 335

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

337 338 339
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
340
            FLOG(FTic counterTimeLevel);
341 342
            const int fackLevel = idxLevel + offsetRealTree;

343 344
            // for each cells
            do{
345
                FLOG(computationCounter.tic());
346
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
347
                FLOG(computationCounter.tac());
348 349 350 351
            } while(octreeIterator.moveRight());

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

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

358 359 360 361 362 363 364 365

    }

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

    /** P2P */
366
    void directPass(const bool p2pEnabled, const bool l2pEnabled){
367
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
368 369 370
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
371 372

        const int heightMinusOne = OctreeHeight - 1;
373
        const FReal boxWidth = tree->getBoxWidth();
374 375 376 377

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
378 379 380
        ContainerClass* neighbors[26];
        FTreeCoordinate offsets[26];
        int neighborPositions[26];
381
        bool hasPeriodicLeaves;
382 383
        // for each leafs
        do{
384 385 386 387 388 389 390 391
            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();
392
                const int counter = tree->getPeriodicLeafsNeighbors( neighbors, neighborPositions, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
393 394 395
                int periodicNeighborsCounter = 0;

                if(hasPeriodicLeaves){
396 397
                    ContainerClass* periodicNeighbors[26];
                    int periodicNeighborPositions[26];
398

399 400
                    for(int idxNeig = 0 ; idxNeig < counter ; ++idxNeig){
                        if( !offsets[idxNeig].equals(0,0,0) ){
401
                            FAssertLF(octreeIterator.getCurrentListTargets() != neighbors[idxNeig]);
402
                            // Put periodic neighbors into other array
BRAMAS Berenger's avatar
BRAMAS Berenger committed
403 404 405
                            FReal*const positionsX = neighbors[idxNeig]->getWPositions()[0];
                            FReal*const positionsY = neighbors[idxNeig]->getWPositions()[1];
                            FReal*const positionsZ = neighbors[idxNeig]->getWPositions()[2];
406

BRAMAS Berenger's avatar
BRAMAS Berenger committed
407
                            for(FSize idxPart = 0; idxPart < neighbors[idxNeig]->getNbParticles() ; ++idxPart){
408 409 410 411
                                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
412

413
                            offsets[periodicNeighborsCounter] = offsets[idxNeig];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
414 415 416
                            periodicNeighbors[periodicNeighborsCounter] = neighbors[idxNeig];
                            periodicNeighborPositions[periodicNeighborsCounter] = neighborPositions[idxNeig];
                            ++periodicNeighborsCounter;
417
                        }
418 419 420 421
                        else{
                            neighbors[idxNeig-periodicNeighborsCounter] = neighbors[idxNeig];
                            neighborPositions[idxNeig-periodicNeighborsCounter] = neighborPositions[idxNeig];
                        }
422
                    }
423

424 425
                    FLOG(computationCounterP2P.tic());
                    kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
426
                                 octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborPositions, periodicNeighborsCounter);
427 428
                    FLOG(computationCounterP2P.tac());

429 430 431 432
                    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];
433

434 435 436 437
                        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());
438 439 440
                        }
                    }
                }
441

442 443
                FLOG(computationCounterP2P.tic());
                kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
444
                             octreeIterator.getCurrentListSrc(), neighbors, neighborPositions, counter - periodicNeighborsCounter);
445 446
                FLOG(computationCounterP2P.tac());
            }
447 448 449
        } while(octreeIterator.moveRight());


450 451 452
        FLOG( FLog::Controller << "\tFinished (@Direct Pass (L2P + P2P) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation L2P : " << computationCounterL2P.cumulated() << " s\n" );
        FLOG( FLog::Controller << "\t\t Computation P2P : " << computationCounterP2P.cumulated() << " s\n" );
453 454 455 456 457 458 459

    }

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

460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
    /** 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
      */
481
    void processPeriodicLevels(){
482
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
483
        FLOG(FTic counterTime);
484

BRAMAS Berenger's avatar
BRAMAS Berenger committed
485 486 487 488 489 490 491 492
        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);
493
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
494 495 496 497 498
            {
                CellClass* virtualChild[8];
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
                    FMemUtils::setall(virtualChild,&upperCells[idxLevel],8);
                    kernels->M2M( &upperCells[idxLevel-1], virtualChild, idxLevel);
499
                }
500
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
501
            CellClass*const downerCells = new CellClass[offsetRealTree];
502

BRAMAS Berenger's avatar
BRAMAS Berenger committed
503 504
            {
                const int idxUpperLevel = 2;
505

506 507
                const CellClass* neighbors[342];
                int neighborPositions[342];
508
                int counter = 0;
509 510 511
                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
512
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
513 514
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
515
                                ++counter;
516
                            }
517 518 519
                        }
                    }
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
520
                // compute M2L
521
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, neighborPositions, counter, idxUpperLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
522
            }
523

BRAMAS Berenger's avatar
BRAMAS Berenger committed
524
            for(int idxUpperLevel = 3 ; idxUpperLevel <= offsetRealTree ; ++idxUpperLevel){
525 526
                const CellClass* neighbors[342];
                int neighborPositions[342];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
527 528 529 530 531
                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){
532 533
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
534
                                ++counter;
535 536 537
                            }
                        }
                    }
538 539
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
540
                // compute M2L
541
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, neighborPositions, counter, idxUpperLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
542
            }
543

BRAMAS Berenger's avatar
BRAMAS Berenger committed
544
            {
545 546
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
547
                for(int idxLevel = 2 ; idxLevel < offsetRealTree-1  ; ++idxLevel){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
548 549 550
                    virtualChild[0] = &downerCells[idxLevel];
                    kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
                }
551
            }
552

553 554 555 556 557 558 559 560
            {
                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
561 562 563 564 565
            // L2L from 0 to level 1
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->L2L( &downerCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
566 567
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
568 569
            delete[] upperCells;
            delete[] downerCells;
570
        }
571

572
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
573 574 575
    }


576 577 578 579
};


#endif // FFMMALGORITHMPERIODIC_HPP