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