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

5 6
#include <array>
#include <algorithm>
7 8

#include "../Utils/FGlobal.hpp"
9
#include "../Utils/FGlobalPeriodic.hpp"
10
#include "../Utils/FAssert.hpp"
11
#include "../Utils/FLog.hpp"
12
#include "../Utils/FAlgorithmTimers.hpp"
13

14
#include "../Utils/FTic.hpp"
15
#include "../Utils/FMemUtils.hpp"
16 17 18 19

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

20
#include "FCoreCommon.hpp"
21 22 23 24 25 26 27

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

36
    using multipole_t       = typename CellClass::multipole_t;
37
    using local_expansion_t = typename CellClass::local_expansion_t;
38
    using symbolic_data_t   = CellClass;
39

40 41
    OctreeClass* const tree;        //< The octree to work on
    KernelClass* kernels;           //< The kernels
42

COULAUD Olivier's avatar
COULAUD Olivier committed
43
    const int OctreeHeight;         //< The height of the octree (real height)
44
    const int nbLevelsAboveRoot;    //< The nb of level the user ask to go above the tree (>= -1)
45
    const int offsetRealTree;       //< nbLevelsAboveRoot GetFakeLevel
46

47
    const int leafLevelSeperationCriteria;
48 49 50 51 52 53

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
54
      * @param inUpperLevel this parameter defines the behavior of the periodicity refer to the main doc
55
      *
56
      */
57
    FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inLeafLevelSeperationCriteria = 1)
58
        : tree(inTree) , kernels(nullptr), OctreeHeight(tree->getHeight()),
59
          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 2), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
60

61 62
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
63
        FAssertLF(leafLevelSeperationCriteria < 3, "Separation criteria should be < 3");
64

65 66
        FAbstractAlgorithm::setNbLevelsInTree(extendedTreeHeight());

67
        FLOG(FLog::Controller << "FFmmAlgorithmPeriodic\n");
68 69 70 71 72 73
    }

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

74 75 76 77
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

78 79 80 81 82 83 84 85



    long long int theoricalRepetition() const {
        if( nbLevelsAboveRoot == -1 ){
            // we know it is 3 (-1;+1)
            return 3;
        }
86
        return 6 * (1 << nbLevelsAboveRoot);
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    }


    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 {
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
        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));
        }
121 122 123 124 125 126 127 128 129
    }

    /** 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
      */
130
    FPoint<FReal> extendedBoxCenter() const {
131 132 133 134 135 136 137 138 139 140 141 142
        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();
143

144 145
            const FReal offset = extendedBoxWidth()/FReal(2.0);
            return FPoint<FReal>( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
146 147
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
        }
    }

    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);
        }
169 170 171 172 173 174 175 176 177 178 179 180 181
    }

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

182 183 184 185 186
    int extendedTreeHeightBoundary() const {
        // The real height
        return OctreeHeight + offsetRealTree + 1;
    }

187
protected:
188 189 190 191
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
192
    void executeCore(const unsigned operationsToProceed) override {
193
        FAssertLF(kernels, "kernels cannot be null");
194

195 196 197 198
        tree->forEachCell([this](CellClass* node){
                node->setLevel(node->getLevel() + offsetRealTree);
            });

199
        if(operationsToProceed & FFmmP2M) bottomPass();
200

201
        if(operationsToProceed & FFmmM2M) upwardPass();
202

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

209
        if(operationsToProceed & FFmmL2L) downardPass();
210

211
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass((operationsToProceed & FFmmP2P),(operationsToProceed & FFmmL2P));
212 213 214 215 216

        tree->forEachCell([this](CellClass* node){
                node->setLevel(node->getLevel() - offsetRealTree);
            });

217 218
    }

219

220 221 222 223 224 225
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

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

        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
237 238
            multipole_t* const leaf_multipole
                = &(octreeIterator.getCurrentCell()->getMultipoleData());
239
            const symbolic_data_t* const leaf_symbolic
240
                = octreeIterator.getCurrentCell();
241
            FLOG(computationCounter.tic());
242 243 244
            kernels->P2M(leaf_multipole,
                         leaf_symbolic,
                         octreeIterator.getCurrentListSrc());
245
            FLOG(computationCounter.tac());
246 247
        } while(octreeIterator.moveRight());

248
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << " s)\n" );
249
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
250 251 252 253 254 255 256 257
    }

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

    /** M2M */
    void upwardPass(){
258
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
259 260
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
261 262 263 264 265 266 267 268 269 270

        // 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 ){
271
            FLOG(FTic counterTimeLevel);
272 273 274 275
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
276 277
                multipole_t* const parent_multipole
                    = &(octreeIterator.getCurrentCell()->getMultipoleData());
278
                const symbolic_data_t* const parent_symbolic
279 280 281 282 283 284 285 286 287
                    = octreeIterator.getCurrentCell();

                CellClass** children = octreeIterator.getCurrentChildren();
                std::array<const multipole_t*, 8> child_multipoles;
                std::transform(children, children+8, child_multipoles.begin(),
                               [](CellClass* c) {
                                   return (c == nullptr ? nullptr
                                           : &(c->getMultipoleData()));
                               });
288
                std::array<const symbolic_data_t*, 8> child_symbolics;
289 290
                std::transform(children, children+8, child_symbolics.begin(),
                               [](CellClass* c) {return c;});
291
                FLOG(computationCounter.tic());
292 293 294 295
                kernels->M2M(parent_multipole,
                             parent_symbolic,
                             child_multipoles.data(),
                             child_symbolics.data());
296
                FLOG(computationCounter.tac());
297 298
            } while(octreeIterator.moveRight());

299
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << (idxLevel + offsetRealTree) << ", " << octreeIterator.getCurrentCell()->getLevel() << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
300 301 302 303 304
            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
        }


305
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << " s)\n" );
306
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
307 308 309
    }

    /////////////////////////////////////////////////////////////////////////////
310
    // Transfer
311 312 313
    /////////////////////////////////////////////////////////////////////////////

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

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

322 323
        const CellClass* neighbors[342];
        int neighborPositions[342];
324 325 326

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
327
            FLOG(FTic counterTimeLevel);
328
            const int fakeLevel = idxLevel + offsetRealTree;
329
            const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
330 331
            // for each cells
            do{
332
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
333 334 335 336 337

                if(counter == 0) {
                    continue;
                }

338
                local_expansion_t* const target_local_exp
339
                    = &(octreeIterator.getCurrentCell()->getLocalExpansionData());
340
                const symbolic_data_t* const target_symbolic
341 342
                    = octreeIterator.getCurrentCell();
                std::array<const multipole_t*, 342> neighbor_multipoles;
343
                std::array<const symbolic_data_t*, 342> neighbor_symbolics;
344 345 346 347 348 349 350 351 352
                std::transform(neighbors, neighbors+counter, neighbor_multipoles.begin(),
                               [](const CellClass* c) {
                                   return (c == nullptr ? nullptr
                                           : &(c->getMultipoleData()));
                               });
                std::transform(neighbors, neighbors+counter, neighbor_symbolics.begin(),
                               [](const CellClass* c) {return c;});


353
                FLOG(computationCounter.tic());
354 355 356 357 358 359 360
                kernels->M2L(
                    target_local_exp,
                    target_symbolic,
                    neighbor_multipoles.data(),
                    neighbor_symbolics.data(),
                    neighborPositions,
                    counter);
361
                FLOG(computationCounter.tac());
362
            } while(octreeIterator.moveRight());
363
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fakeLevel << ", " << octreeIterator.getCurrentCell()->getLevel() << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
364 365
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
366

367
            FLOG(computationCounter.tic());
368
            kernels->finishedLevelM2L(fakeLevel);
369
            FLOG(computationCounter.tac());
370
        }
371
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << " s)\n" );
372
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
373
    }
374

375 376 377 378 379 380
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
381
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
382 383
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
384 385 386

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

388 389 390
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
391
            FLOG(FTic counterTimeLevel);
392

393 394
            // for each cells
            do{
395
                local_expansion_t* const parent_local_exp
396
                    = &(octreeIterator.getCurrentCell()->getLocalExpansionData());
397
                const symbolic_data_t* const parent_symbolic
398 399
                    = octreeIterator.getCurrentCell();
                CellClass** children = octreeIterator.getCurrentChildren();
400
                std::array<local_expansion_t*, 8> child_local_expansions;
401 402 403 404
                std::transform(children, children+8, child_local_expansions.begin(),
                               [](CellClass* c) {return (c == nullptr ? nullptr
                                                         : &(c->getLocalExpansionData()));
                               });
405
                std::array<symbolic_data_t*, 8> child_symbolics;
406 407
                std::transform(children, children+8, child_symbolics.begin(),
                               [](CellClass* c) {return c;});
408
                FLOG(computationCounter.tic());
409 410 411 412 413 414
                kernels->L2L(
                    parent_local_exp,
                    parent_symbolic,
                    child_local_expansions.data(),
                    child_symbolics.data()
                    );
415
                FLOG(computationCounter.tac());
416 417
            } while(octreeIterator.moveRight());

418
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << (idxLevel + offsetRealTree) << ", " << octreeIterator.getCurrentCell()->getLevel() << ") = "  << counterTimeLevel.tacAndElapsed() << " s\n" );
419 420
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
421 422
        }

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

426 427 428 429 430 431 432 433

    }

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

    /** P2P */
434
    void directPass(const bool p2pEnabled, const bool l2pEnabled){
435
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
436 437 438
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
439 440

        const int heightMinusOne = OctreeHeight - 1;
441
        const FReal boxWidth = tree->getBoxWidth();
442 443 444 445

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
446 447 448
        ContainerClass* neighbors[26];
        FTreeCoordinate offsets[26];
        int neighborPositions[26];
449
        bool hasPeriodicLeaves;
450 451
        // for each leafs
        do{
452 453
            if(l2pEnabled){
                FLOG(computationCounterL2P.tic());
454 455 456 457 458
                kernels->L2P(
                    &(octreeIterator.getCurrentCell()->getLocalExpansionData()),
                    octreeIterator.getCurrentCell(),
                    octreeIterator.getCurrentListTargets()
                    );
459 460 461 462 463
                FLOG(computationCounterL2P.tac());
            }
            if(p2pEnabled){
                // need the current particles and neighbors particles
                const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
464 465 466
                const int counter = tree->getPeriodicLeafsNeighbors( neighbors, neighborPositions, offsets,
                                                                     hasPeriodicLeaves, centerOfLeaf,
                                                                     heightMinusOne, AllDirs);
467 468 469
                int periodicNeighborsCounter = 0;

                if(hasPeriodicLeaves){
470 471
                    ContainerClass* periodicNeighbors[26];
                    int periodicNeighborPositions[26];
472

473 474
                    for(int idxNeig = 0 ; idxNeig < counter ; ++idxNeig){
                        if( !offsets[idxNeig].equals(0,0,0) ){
475
                            FAssertLF(octreeIterator.getCurrentListTargets() != neighbors[idxNeig]);
476
                            // Put periodic neighbors into other array
477 478 479
                            FReal*const positionsX = neighbors[idxNeig]->getPositions()[0];
                            FReal*const positionsY = neighbors[idxNeig]->getPositions()[1];
                            FReal*const positionsZ = neighbors[idxNeig]->getPositions()[2];
480

BRAMAS Berenger's avatar
BRAMAS Berenger committed
481
                            for(FSize idxPart = 0; idxPart < neighbors[idxNeig]->getNbParticles() ; ++idxPart){
482 483 484 485
                                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
486

487
                            offsets[periodicNeighborsCounter] = offsets[idxNeig];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
488 489 490
                            periodicNeighbors[periodicNeighborsCounter] = neighbors[idxNeig];
                            periodicNeighborPositions[periodicNeighborsCounter] = neighborPositions[idxNeig];
                            ++periodicNeighborsCounter;
491
                        }
492 493 494 495
                        else{
                            neighbors[idxNeig-periodicNeighborsCounter] = neighbors[idxNeig];
                            neighborPositions[idxNeig-periodicNeighborsCounter] = neighborPositions[idxNeig];
                        }
496
                    }
497

498 499
                    FLOG(computationCounterP2P.tic());
                    kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
500
                                 octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborPositions, periodicNeighborsCounter);
501 502
                    FLOG(computationCounterP2P.tac());

503
                    for(int idxNeig = 0 ; idxNeig < periodicNeighborsCounter ; ++idxNeig){
504 505 506
                        FReal*const positionsX = periodicNeighbors[idxNeig]->getPositions()[0];
                        FReal*const positionsY = periodicNeighbors[idxNeig]->getPositions()[1];
                        FReal*const positionsZ = periodicNeighbors[idxNeig]->getPositions()[2];
507

508 509 510 511
                        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());
512 513 514
                        }
                    }
                }
515

516 517
                FLOG(computationCounterP2P.tic());
                kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
518
                             octreeIterator.getCurrentListSrc(), neighbors, neighborPositions, counter - periodicNeighborsCounter);
519 520
                FLOG(computationCounterP2P.tac());
            }
521 522 523
        } while(octreeIterator.moveRight());


524
        FLOG( FLog::Controller << "\tFinished (@Direct Pass (L2P + P2P) = "  << counterTime.tacAndElapsed() << " s)\n" );
525 526
        FLOG( FLog::Controller << "\t\t Computation L2P : " << computationCounterL2P.cumulated() << " s\n" );
        FLOG( FLog::Controller << "\t\t Computation P2P : " << computationCounterP2P.cumulated() << " s\n" );
527 528 529 530 531 532 533

    }

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

534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
    /** 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
      */
555
    void processPeriodicLevels(){
556
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
557
        FLOG(FTic counterTime);
558

BRAMAS Berenger's avatar
BRAMAS Berenger committed
559 560 561 562
        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];
563 564 565 566
            for(int i = 0; i < offsetRealTree; ++i) {
                upperCells[i].setLevel(i+1); // TODO: check that the level is correct (see other todos under)
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
567 568 569
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
570 571 572

                multipole_t* const real_tree_root_multipole
                    = &((&upperCells[offsetRealTree-1])->getMultipoleData());
573
                const symbolic_data_t* const real_tree_root_symbolic
574 575 576 577 578 579 580 581 582 583
                    = &(upperCells[offsetRealTree-1]);

                CellClass** children = octreeIterator.getCurrentBox();

                std::array<const multipole_t*, 8> level_1_multipoles;
                std::transform(children, children+8, level_1_multipoles.begin(),
                               [](CellClass* c) {
                                   return (c == nullptr ? nullptr
                                           : &(c->getMultipoleData()));
                               });
584
                std::array<const symbolic_data_t*, 8> level_1_symbolics;
585 586 587 588 589 590 591 592 593
                std::transform(children, children+8, level_1_symbolics.begin(),
                               [](CellClass* c) {return c;});
                // TODO check that parent symbolic has the right level
                kernels->M2M(
                    real_tree_root_multipole,
                    real_tree_root_symbolic,
                    level_1_multipoles.data(),
                    level_1_symbolics.data()
                    );
594
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
595
            {
596
                CellClass* children[8] = {};
BRAMAS Berenger's avatar
BRAMAS Berenger committed
597
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
598 599 600 601 602
                    // Copy virtual cell at given level into virtual children
                    FMemUtils::setall(children,&upperCells[idxLevel],8);

                    multipole_t* const virtual_parent_multipole
                        = &((&upperCells[idxLevel-1])->getMultipoleData());
603
                    const symbolic_data_t* const virtual_parent_symbolic
604 605 606 607 608
                        = &(upperCells[idxLevel-1]);

                    std::array<const multipole_t*, 8> virtual_child_multipoles;
                    std::transform(children, children+8, virtual_child_multipoles.begin(),
                                   [](CellClass* c) {return &(c->getMultipoleData());});
609
                    std::array<const symbolic_data_t*, 8> virtual_child_symbolics;
610 611 612 613 614 615 616 617 618
                    std::transform(children, children+8, virtual_child_symbolics.begin(),
                                   [](CellClass* c) {return c;});

                    kernels->M2M(
                        virtual_parent_multipole,
                        virtual_parent_symbolic,
                        virtual_child_multipoles.data(),
                        virtual_child_symbolics.data()
                        );
619
                }
620
            }
621

622 623 624 625 626
            CellClass*const downerCells = new CellClass[offsetRealTree];
            for(int i = 0; i < offsetRealTree; ++i) {
                downerCells[i].setLevel(i+1); // TODO: check that the level is correct (see other todos under)
            }
            // Build a virtual environment for the M2L at the topmost level
BRAMAS Berenger's avatar
BRAMAS Berenger committed
627 628
            {
                const int idxUpperLevel = 2;
629

630 631
                const CellClass* neighbors[342];
                int neighborPositions[342];
632
                int counter = 0;
633 634 635

                // The cells are all the same, duplicate the central cell and
                // create fake neighbour positions
636 637 638
                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
639
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
640 641
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
642
                                ++counter;
643
                            }
644 645 646
                        }
                    }
                }
647

648
                local_expansion_t* const target_local_exp
649
                    = &(downerCells[idxUpperLevel-1].getLocalExpansionData());
650
                const symbolic_data_t* const target_symbolic
651 652
                    = &downerCells[idxUpperLevel-1];
                std::array<const multipole_t*, 342> neighbor_multipoles;
653
                std::array<const symbolic_data_t*, 342> neighbor_symbolics;
654 655 656 657 658
                std::transform(neighbors, neighbors+counter, neighbor_multipoles.begin(),
                               [](const CellClass* c) {return &(c->getMultipoleData());});
                std::transform(neighbors, neighbors+counter, neighbor_symbolics.begin(),
                               [](const CellClass* c) {return c;});

BRAMAS Berenger's avatar
BRAMAS Berenger committed
659
                // compute M2L
660 661 662 663 664 665 666
                kernels->M2L(
                    target_local_exp,
                    target_symbolic,
                    neighbor_multipoles.data(),
                    neighbor_symbolics.data(),
                    neighborPositions,
                    counter);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
667
            }
668 669
            // note: the triple loop bounds are not the same than for the previous piece of code
            // which handles the topmost virtual cell
BRAMAS Berenger's avatar
BRAMAS Berenger committed
670
            for(int idxUpperLevel = 3 ; idxUpperLevel <= offsetRealTree ; ++idxUpperLevel){
671 672
                const CellClass* neighbors[342];
                int neighborPositions[342];
BRAMAS Berenger's avatar
BRAMAS Berenger committed
673 674 675 676 677
                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){
678 679
                                neighbors[counter] = &upperCells[idxUpperLevel-1];
                                neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
680
                                ++counter;
681 682 683
                            }
                        }
                    }
684 685
                }

686
                local_expansion_t* const target_local_exp
687
                    = &(downerCells[idxUpperLevel-1].getLocalExpansionData());
688
                const symbolic_data_t* const target_symbolic
689 690
                    = &downerCells[idxUpperLevel-1];
                std::array<const multipole_t*, 342> neighbor_multipoles;
691
                std::array<const symbolic_data_t*, 342> neighbor_symbolics;
692 693 694 695 696 697 698 699
                std::transform(neighbors, neighbors+counter, neighbor_multipoles.begin(),
                               [](const CellClass* c) {
                                   return (c == nullptr ? nullptr
                                           : &(c->getMultipoleData()));
                               });
                std::transform(neighbors, neighbors+counter, neighbor_symbolics.begin(),
                               [](const CellClass* c) {return c;});

BRAMAS Berenger's avatar
BRAMAS Berenger committed
700
                // compute M2L
701 702 703 704 705 706 707
                kernels->M2L(
                    target_local_exp,
                    target_symbolic,
                    neighbor_multipoles.data(),
                    neighbor_symbolics.data(),
                    neighborPositions,
                    counter);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
708
            }
709

710 711

            // Run the L2L for all but the lowest virtual levels
BRAMAS Berenger's avatar
BRAMAS Berenger committed
712
            {
713 714
                std::array<local_expansion_t*, 8> virtual_child_local_exps = {};
                std::array<symbolic_data_t*, 8> virtual_child_symbolics = {};
715
                for(int idxLevel = 2 ; idxLevel < offsetRealTree-1  ; ++idxLevel){
716
                    local_expansion_t* const virtual_parent_local_exp
717
                        = &(downerCells[idxLevel-1].getLocalExpansionData());
718
                    const symbolic_data_t* const virtual_parent_symbolic
719 720 721 722 723 724 725 726 727 728 729 730 731
                        = &downerCells[idxLevel-1];

                    virtual_child_local_exps[0]
                        = &(downerCells[idxLevel].getLocalExpansionData());
                    virtual_child_symbolics[0]
                        = &(downerCells[idxLevel]);

                    kernels->L2L(
                        virtual_parent_local_exp,
                        virtual_parent_symbolic,
                        virtual_child_local_exps.data(),
                        virtual_child_symbolics.data()
                        );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
732
                }
733
            }
734
            // Run the L2L for the lowest virtual level
735
            {
736 737
                std::array<local_expansion_t*, 8> virtual_child_local_exps = {};
                std::array<symbolic_data_t*, 8> virtual_child_symbolics = {};
738
                const int idxLevel = offsetRealTree-1;
739
                local_expansion_t* const virtual_parent_local_exp
740
                    = &(downerCells[idxLevel-1].getLocalExpansionData());
741
                const symbolic_data_t* const virtual_parent_symbolic
742 743 744 745 746 747 748 749 750
                    = &downerCells[idxLevel-1];
                virtual_child_local_exps[7] = &(downerCells[idxLevel].getLocalExpansionData());
                virtual_child_symbolics[7] = &downerCells[idxLevel];
                kernels->L2L(
                    virtual_parent_local_exp,
                    virtual_parent_symbolic,
                    virtual_child_local_exps.data(),
                    virtual_child_symbolics.data()
                    );
751 752
            }

753
            // Run L2L from the lowest vitual level to the highest real tree level
BRAMAS Berenger's avatar
BRAMAS Berenger committed
754 755 756
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
757

758
                local_expansion_t* const virtual_parent_local_exp
759
                    = &(downerCells[offsetRealTree-1].getLocalExpansionData());
760
                const symbolic_data_t* const virtual_parent_symbolic
761 762 763
                    = &downerCells[offsetRealTree-1];

                CellClass** children = octreeIterator.getCurrentBox();
764
                std::array<local_expansion_t*, 8> child_local_expansions;
765 766 767 768
                std::transform(children, children+8, child_local_expansions.begin(),
                               [](CellClass* c) {return (c == nullptr ? nullptr
                                                         : &(c->getLocalExpansionData()));
                               });
769
                std::array<symbolic_data_t*, 8> child_symbolics;
770 771 772 773 774 775 776 777 778
                std::transform(children, children+8, child_symbolics.begin(),
                               [](CellClass* c) {return c;});

                kernels->L2L(
                    virtual_parent_local_exp,
                    virtual_parent_symbolic,
                    child_local_expansions.data(),
                    child_symbolics.data()
                    );
779 780
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
781 782
            delete[] upperCells;
            delete[] downerCells;
783
        }
784

785
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << " s)\n" );
786 787 788
    }


789 790 791 792
};


#endif // FFMMALGORITHMPERIODIC_HPP