Une MAJ de sécurité est nécessaire sur notre version actuelle. Elle sera effectuée lundi 02/08 entre 12h30 et 13h. L'interruption de service devrait durer quelques minutes (probablement moins de 5 minutes).

FFmmAlgorithmPeriodic.hpp 22.1 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"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
22
#include "../Utils/FAssert.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
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 OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
BRAMAS Berenger's avatar
BRAMAS Berenger committed
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 55 56 57 58 59

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

BRAMAS Berenger's avatar
BRAMAS Berenger committed
67 68
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
69

BRAMAS Berenger's avatar
BRAMAS Berenger committed
70
        FLOG(FLog::Controller << "FFmmAlgorithmPeriodic\n");
71 72 73 74 75 76
    }

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

77 78 79 80
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

81 82 83 84
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
85
    void execute(const unsigned operationsToProceed = FFmmNearAndFarFields){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
86
        FAssertLF(kernels, "kernels cannot be null");
87

88
        if(operationsToProceed & FFmmP2M) bottomPass();
89

90
        if(operationsToProceed & FFmmM2M) upwardPass();
91

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

98
        if(operationsToProceed & FFmmL2L) downardPass();
99

100
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass((operationsToProceed & FFmmP2P),(operationsToProceed & FFmmL2P));
101 102
    }

103

104 105 106 107 108 109
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
110
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
111 112
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
113 114 115 116 117 118 119 120

        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
BRAMAS Berenger's avatar
BRAMAS Berenger committed
121
            FLOG(computationCounter.tic());
122
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
BRAMAS Berenger's avatar
BRAMAS Berenger committed
123
            FLOG(computationCounter.tac());
124 125
        } while(octreeIterator.moveRight());

BRAMAS Berenger's avatar
BRAMAS Berenger committed
126 127
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
128 129 130 131 132 133 134 135
    }

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

    /** M2M */
    void upwardPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
136
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
137 138
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
139 140 141 142 143 144 145 146 147 148

        // 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 ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
149
            FLOG(FTic counterTimeLevel);
150
            const int fackLevel = idxLevel + offsetRealTree;
151 152 153 154
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
BRAMAS Berenger's avatar
BRAMAS Berenger committed
155
                FLOG(computationCounter.tic());
156
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
157
                FLOG(computationCounter.tac());
158 159 160 161
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
162
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
163 164 165
        }


BRAMAS Berenger's avatar
BRAMAS Berenger committed
166 167
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
168 169 170
    }

    /////////////////////////////////////////////////////////////////////////////
171
    // Transfer
172 173 174
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
175
    void transferPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
176
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
177 178
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
179 180 181 182

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

183
        const CellClass* neighbors[343];
184 185 186

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
187
            FLOG(FTic counterTimeLevel);
188
            const int fackLevel = idxLevel + offsetRealTree;
189 190
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
191
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
192
                FLOG(computationCounter.tic());
193
                if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
194
                FLOG(computationCounter.tac());
195 196 197
            } while(octreeIterator.moveRight());
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
198

BRAMAS Berenger's avatar
BRAMAS Berenger committed
199
            FLOG(computationCounter.tic());
200
            kernels->finishedLevelM2L(fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
201
            FLOG(computationCounter.tac());
BRAMAS Berenger's avatar
BRAMAS Berenger committed
202
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
203
        }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
204 205
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
206
    }
207

208 209 210 211 212 213
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
BRAMAS Berenger's avatar
BRAMAS Berenger committed
214
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
215 216
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
217 218 219

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

221 222 223
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
224
            FLOG(FTic counterTimeLevel);
225 226
            const int fackLevel = idxLevel + offsetRealTree;

227 228
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
229
                FLOG(computationCounter.tic());
230
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
231
                FLOG(computationCounter.tac());
232 233 234 235
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
236
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
237 238
        }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
239 240
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (L2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
241

242 243 244 245 246 247 248 249

    }

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

    /** P2P */
250
    void directPass(const bool p2pEnabled, const bool l2pEnabled){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
251
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
252 253 254
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
255 256

        const int heightMinusOne = OctreeHeight - 1;
257
        const FReal boxWidth = tree->getBoxWidth();
258 259 260 261

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
berenger-bramas's avatar
berenger-bramas committed
262
        ContainerClass* neighbors[27];
263
        FTreeCoordinate offsets[27];
264
        bool hasPeriodicLeaves;
265 266
        // for each leafs
        do{
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
            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();
                const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
                int periodicNeighborsCounter = 0;

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

                    for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
                        if( neighbors[idxNeig] && !offsets[idxNeig].equals(0,0,0) ){
                            // Put periodic neighbors into other array
                            periodicNeighbors[idxNeig] = neighbors[idxNeig];
                            neighbors[idxNeig] = nullptr;
                            ++periodicNeighborsCounter;

                            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());
                            }
298 299
                        }
                    }
300

301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316
                    FLOG(computationCounterP2P.tic());
                    kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                                 octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborsCounter);
                    FLOG(computationCounterP2P.tac());

                    for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
                        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());
                            }
317 318 319
                        }
                    }
                }
320

321 322 323 324 325
                FLOG(computationCounterP2P.tic());
                kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                             octreeIterator.getCurrentListSrc(), neighbors, counter - periodicNeighborsCounter);
                FLOG(computationCounterP2P.tac());
            }
326 327 328
        } while(octreeIterator.moveRight());


BRAMAS Berenger's avatar
BRAMAS Berenger committed
329 330 331
        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" );
332 333 334 335 336 337 338

    }

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

339 340 341 342 343 344 345 346 347 348
    /** 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);
    }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
349

350
    long long int theoricalRepetition() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
351 352 353 354 355 356 357 358
        if( nbLevelsAboveRoot == -1 ){
            // we know it is 3 (-1;+1)
            return 3;
        }
        // Else we find the repetition in one dir and double it
        const long long int oneDirectionRepetition = (1<<(nbLevelsAboveRoot+2)); // 2^nbLevelsAboveRoot in each dim
        const long long int fullRepetition = 2 * oneDirectionRepetition;
        return fullRepetition;
359 360
    }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
361

362
    void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
363 364 365 366 367 368 369 370 371 372 373
        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);
        }
374 375 376
    }


377
    FReal extendedBoxWidth() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
378 379 380
        // The simulation box is repeated is repeated 4 times if nbLevelsAboveRoot==-1
        // And then it doubles by two
        return tree->getBoxWidth() * FReal(1<<(nbLevelsAboveRoot+3));
381 382
    }

383 384 385 386 387 388 389 390
    /** 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 {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
391
        const FReal originalBoxWidth     = tree->getBoxWidth();
392
        const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
393
        const FPoint originalBoxCenter   = tree->getBoxCenter();
394

BRAMAS Berenger's avatar
BRAMAS Berenger committed
395
        const FReal offset = extendedBoxWidth()/FReal(2.0);
396 397 398
        return FPoint( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
399 400 401 402 403 404 405 406
    }

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

412 413 414 415 416 417 418 419 420 421 422
    /** 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
      */
423
    void processPeriodicLevels(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
424
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
425
        FLOG(FTic counterTime);
426

BRAMAS Berenger's avatar
BRAMAS Berenger committed
427 428 429 430 431 432 433 434
        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);
435
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
436 437 438 439 440
            {
                CellClass* virtualChild[8];
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
                    FMemUtils::setall(virtualChild,&upperCells[idxLevel],8);
                    kernels->M2M( &upperCells[idxLevel-1], virtualChild, idxLevel);
441
                }
442
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
443
            CellClass*const downerCells = new CellClass[offsetRealTree];
444

BRAMAS Berenger's avatar
BRAMAS Berenger committed
445 446
            {
                const int idxUpperLevel = 2;
447 448 449 450

                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                int counter = 0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
451 452 453
                for(int idxX = -2 ; idxX <= 1 ; ++idxX){
                    for(int idxY = -2 ; idxY <= 1 ; ++idxY){
                        for(int idxZ = -2 ; idxZ <= 1 ; ++idxZ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
454 455 456
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
                                neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
                                ++counter;
457
                            }
458 459 460
                        }
                    }
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
461 462 463
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
464

BRAMAS Berenger's avatar
BRAMAS Berenger committed
465 466 467 468 469 470 471 472 473 474
            for(int idxUpperLevel = 3 ; idxUpperLevel <= offsetRealTree ; ++idxUpperLevel){
                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                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){
                                neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
                                ++counter;
475 476 477
                            }
                        }
                    }
478 479
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
480 481 482
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
483

BRAMAS Berenger's avatar
BRAMAS Berenger committed
484
            {
485 486
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
487 488 489 490
                for(int idxLevel = 2 ; idxLevel <= offsetRealTree-1  ; ++idxLevel){
                    virtualChild[0] = &downerCells[idxLevel];
                    kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
                }
491
            }
492

BRAMAS Berenger's avatar
BRAMAS Berenger committed
493 494 495 496 497
            // L2L from 0 to level 1
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->L2L( &downerCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
498 499
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
500 501
            delete[] upperCells;
            delete[] downerCells;
502
        }
503

BRAMAS Berenger's avatar
BRAMAS Berenger committed
504
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
505 506 507
    }


508 509 510 511
};


#endif // FFMMALGORITHMPERIODIC_HPP