FFmmAlgorithmPeriodic.hpp 25.5 KB
Newer Older
1 2
#ifndef FFMMALGORITHMPERIODIC_HPP
#define FFMMALGORITHMPERIODIC_HPP
3

4 5

#include "../Utils/FGlobal.hpp"
6
#include "../Utils/FGlobalPeriodic.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
7
#include "../Utils/FAssert.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
8
#include "../Utils/FLog.hpp"
9

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

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

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

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

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

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

39 40 41 42 43 44

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

BRAMAS Berenger's avatar
BRAMAS Berenger committed
52 53
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
54

BRAMAS Berenger's avatar
BRAMAS Berenger committed
55
        FLOG(FLog::Controller << "FFmmAlgorithmPeriodic\n");
56 57 58 59 60 61
    }

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

62 63 64 65
    void setKernel(KernelClass*const inKernel){
        kernels = inKernel;
    }

66 67 68 69
    /**
      * To execute the fmm algorithm
      * Call this function to run the complete algorithm
      */
70
    void execute(const unsigned operationsToProceed = FFmmNearAndFarFields){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
71
        FAssertLF(kernels, "kernels cannot be null");
72

73
        if(operationsToProceed & FFmmP2M) bottomPass();
74

75
        if(operationsToProceed & FFmmM2M) upwardPass();
76

77 78 79 80 81
        if(operationsToProceed & FFmmM2L){
            transferPass();
            // before downward pass we have to perform the periodicity
            processPeriodicLevels();
        }
82

83
        if(operationsToProceed & FFmmL2L) downardPass();
84

BRAMAS Berenger's avatar
BRAMAS Berenger committed
85
        if((operationsToProceed & FFmmP2P) || (operationsToProceed & FFmmL2P)) directPass();
86 87
    }

88

89 90 91 92 93 94
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
95
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
96 97
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
98 99 100 101 102 103 104 105

        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
106
            FLOG(computationCounter.tic());
107
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
BRAMAS Berenger's avatar
BRAMAS Berenger committed
108
            FLOG(computationCounter.tac());
109 110
        } while(octreeIterator.moveRight());

BRAMAS Berenger's avatar
BRAMAS Berenger committed
111 112
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
113 114 115 116 117 118 119 120
    }

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

    /** M2M */
    void upwardPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
121
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
122 123
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
124 125 126 127 128 129 130 131 132 133

        // 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
134
            FLOG(FTic counterTimeLevel);
135
            const int fackLevel = idxLevel + offsetRealTree;
136 137 138 139
            // 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
140
                FLOG(computationCounter.tic());
141
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
142
                FLOG(computationCounter.tac());
143 144 145 146
            } while(octreeIterator.moveRight());

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


BRAMAS Berenger's avatar
BRAMAS Berenger committed
151 152
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
153 154 155
    }

    /////////////////////////////////////////////////////////////////////////////
156
    // Transfer
157 158 159
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
160
    void transferPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
161
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
162 163
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
164 165 166 167

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

168
        const CellClass* neighbors[343];
169 170 171

        // for each levels
        for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
172
            FLOG(FTic counterTimeLevel);
173
            const int fackLevel = idxLevel + offsetRealTree;
174 175
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
176
                const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
177
                FLOG(computationCounter.tic());
178
                if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
179
                FLOG(computationCounter.tac());
180 181 182
            } while(octreeIterator.moveRight());
            avoidGotoLeftIterator.moveDown();
            octreeIterator = avoidGotoLeftIterator;
183

BRAMAS Berenger's avatar
BRAMAS Berenger committed
184
            FLOG(computationCounter.tic());
185
            kernels->finishedLevelM2L(fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
186
            FLOG(computationCounter.tac());
BRAMAS Berenger's avatar
BRAMAS Berenger committed
187
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
188
        }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
189 190
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
191
    }
192

193 194 195 196 197 198
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
BRAMAS Berenger's avatar
BRAMAS Berenger committed
199
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
200 201
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
202 203 204

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

206 207 208
        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
209
            FLOG(FTic counterTimeLevel);
210 211
            const int fackLevel = idxLevel + offsetRealTree;

212 213
            // for each cells
            do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
214
                FLOG(computationCounter.tic());
215
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
216
                FLOG(computationCounter.tac());
217 218 219 220
            } while(octreeIterator.moveRight());

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

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

227 228 229 230 231 232 233 234 235

    }

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

    /** P2P */
    void directPass(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
236
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
237 238 239
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
240 241

        const int heightMinusOne = OctreeHeight - 1;
242
        const FReal boxWidth = tree->getBoxWidth();
243 244 245 246

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
berenger-bramas's avatar
berenger-bramas committed
247
        ContainerClass* neighbors[27];
248
        FTreeCoordinate offsets[27];
249
        bool hasPeriodicLeaves;
250 251
        // for each leafs
        do{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
252
            FLOG(computationCounterL2P.tic());
253
            kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
BRAMAS Berenger's avatar
BRAMAS Berenger committed
254
            FLOG(computationCounterL2P.tac());
255

256
            // need the current particles and neighbors particles
257
            const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
258
            const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
259 260 261 262 263
            int periodicNeighborsCounter = 0;

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

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
266 267 268
                    if( neighbors[idxNeig] && !offsets[idxNeig].equals(0,0,0) ){
                        // Put periodic neighbors into other array
                        periodicNeighbors[idxNeig] = neighbors[idxNeig];
269
                        neighbors[idxNeig] = nullptr;
270
                        ++periodicNeighborsCounter;
271 272 273 274 275 276 277 278 279

                        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());
280 281 282 283
                        }
                    }
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
284
                FLOG(computationCounterP2P.tic());
285 286
                kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                             octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborsCounter);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
287
                FLOG(computationCounterP2P.tac());
288 289

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
290
                    if( periodicNeighbors[idxNeig] ){
291 292 293 294 295 296 297 298
                        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());
299 300 301 302 303
                        }
                    }
                }
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
304
            FLOG(computationCounterP2P.tic());
305 306
            kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                         octreeIterator.getCurrentListSrc(), neighbors, counter - periodicNeighborsCounter);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
307
            FLOG(computationCounterP2P.tac());
308 309


310 311 312
        } while(octreeIterator.moveRight());


BRAMAS Berenger's avatar
BRAMAS Berenger committed
313 314 315
        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" );
316 317 318 319 320 321 322

    }

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

323

324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
    /** This function process several M2M from level nbLevelsAboveRoot to level 0
      * and give the final result
      * @param result the cell at the last M2M
      * @param root the starting cell
      * @param startX the beginning of the index in x [0;endX]
      * @param endX the end of the index in x [startX;1]
      * @param startY the beginning of the index in y [0;endY]
      * @param endY the end of the index in y [startY;1]
      * @param startZ the beginning of the index in z [0;endZ]
      * @param endZ the end of the index in z [startZ;1]
      */
    void processTopM2MInIntervals( CellClass*const result, const CellClass& root, const int startX,
                              const int endX, const int startY, const int endY, const int startZ,
                              const int endZ){
        // allocate array
        CellClass*const cellsAtLevel = new CellClass[nbLevelsAboveRoot+2];
        // process by using other function
        processM2MInIntervals(cellsAtLevel,root,startX,endX,startY,endY,startZ,endZ);
        // copy result
        *result = cellsAtLevel[0];
        delete[] cellsAtLevel;
    }

    /** This function process several M2M from level nbLevelsAboveRoot to level 0
      * @param cellsAtLevel the intermediate results
      * @param root the starting cell
      * @param startX the beginning of the index in x [0;endX]
      * @param endX the end of the index in x [startX;1]
      * @param startY the beginning of the index in y [0;endY]
      * @param endY the end of the index in y [startY;1]
      * @param startZ the beginning of the index in z [0;endZ]
      * @param endZ the end of the index in z [startZ;1]
      */
    void  processM2MInIntervals( CellClass cellsAtLevel[], const CellClass& root, const int startX,
                              const int endX, const int startY, const int endY, const int startZ,
                              const int endZ){
        // start from the initial cell
        cellsAtLevel[nbLevelsAboveRoot+1] = root;
        // to create virtual children
        CellClass* virtualChild[8];
        // for all levels
        for(int idxLevel = nbLevelsAboveRoot ; idxLevel >= 0  ; --idxLevel){
            // reset children
            memset(virtualChild, 0, sizeof(CellClass*)*8);
            // fill the vector with previous result
            for(int idxX = startX ; idxX <= endX ; ++idxX){
                for(int idxY = startY ; idxY <= endY ; ++idxY){
                    for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
                        virtualChild[childIndex(idxX,idxY,idxZ)] = &cellsAtLevel[idxLevel+1];
                    }
                }
            }
            // compute the M2M
            kernels->M2M( &cellsAtLevel[idxLevel], virtualChild, idxLevel + 2);
        }
    }

    /** Fill an interactions neighbors with some intervals
      * @param neighbors the vector to fill
      * @param source the source cell to fill the vector
      * @param startX the beginning of the index in x [-3;0]
      * @param endX the end of the index in x  [0;3]
      * @param startY the beginning of the index in y [-3;0]
      * @param endY the end of the index in y [0;3]
      * @param startZ the beginning of the index in z [-3;0]
      * @param endZ the end of the index in z [0;3]
      * @return the number of position filled
      */
    int  fillM2LVectorFromIntervals(const CellClass* neighbors[343], const CellClass& source,
                     const int startX, const int endX, const int startY, const int endY,
                     const int startZ, const int endZ){
        int counter = 0;
        // for all x in interval
        for(int idxX = startX ; idxX <= endX ; ++idxX){
            // for all y in interval
            for(int idxY = startY ; idxY <= endY ; ++idxY){
                // for all z in interval
                for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
                    // do not fill close neigbors
BRAMAS Berenger's avatar
BRAMAS Berenger committed
403
                    if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1 ){
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
                        neighbors[neighIndex(idxX,idxY,idxZ)] = &source;
                        ++counter;
                    }
                }
            }
        }
        // return the number of position filled
        return counter;
    }

    /** Get the index of a child (for the M2M and the L2L)
      * @param x the x position in the children  (from 0 to +1)
      * @param y the y position in the children  (from 0 to +1)
      * @param z the z position in the children  (from 0 to +1)
      * @return the index (from 0 to 7)
      */
    int childIndex(const int x, const int y, const int z) const {
        return (x<<2) | (y<<1) | z;
    }

    /** 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
434

435
    long long int theoricalRepetition() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
436 437 438 439 440 441 442 443
        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;
444 445
    }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
446

447
    void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
448 449 450 451 452 453 454 455 456 457 458
        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);
        }
459 460 461
    }


462
    FReal extendedBoxWidth() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
463 464 465
        // 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));
466 467
    }

468 469 470 471 472 473 474 475
    /** 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
476
        const FReal originalBoxWidth     = tree->getBoxWidth();
477
        const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
478
        const FPoint originalBoxCenter   = tree->getBoxCenter();
479

BRAMAS Berenger's avatar
BRAMAS Berenger committed
480
        const FReal offset = extendedBoxWidth()/FReal(2.0);
481 482 483
        return FPoint( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
484 485 486 487 488 489 490 491
    }

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

497 498 499 500 501 502 503 504 505 506 507
    /** 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
      */
508
    void processPeriodicLevels(){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
509
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
BRAMAS Berenger's avatar
BRAMAS Berenger committed
510
        FLOG(FTic counterTime);
511

BRAMAS Berenger's avatar
BRAMAS Berenger committed
512 513 514 515 516 517 518 519
        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);
520
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
521 522 523 524 525
            {
                CellClass* virtualChild[8];
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
                    FMemUtils::setall(virtualChild,&upperCells[idxLevel],8);
                    kernels->M2M( &upperCells[idxLevel-1], virtualChild, idxLevel);
526
                }
527
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
528
            CellClass*const downerCells = new CellClass[offsetRealTree];
529

BRAMAS Berenger's avatar
BRAMAS Berenger committed
530 531
            {
                const int idxUpperLevel = 2;
532 533 534 535

                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                int counter = 0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
536 537 538
                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
539 540 541
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
                                neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
                                ++counter;
542
                            }
543 544 545
                        }
                    }
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
546 547 548
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
549

BRAMAS Berenger's avatar
BRAMAS Berenger committed
550 551 552 553 554 555 556 557 558 559
            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;
560 561 562
                            }
                        }
                    }
563 564
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
565 566 567
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
568

BRAMAS Berenger's avatar
BRAMAS Berenger committed
569
            {
570 571
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
572 573 574 575
                for(int idxLevel = 2 ; idxLevel <= offsetRealTree-1  ; ++idxLevel){
                    virtualChild[0] = &downerCells[idxLevel];
                    kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
                }
576
            }
577

BRAMAS Berenger's avatar
BRAMAS Berenger committed
578 579 580 581 582
            // L2L from 0 to level 1
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->L2L( &downerCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
583 584
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
585 586
            delete[] upperCells;
            delete[] downerCells;
587
        }
588

BRAMAS Berenger's avatar
BRAMAS Berenger committed
589
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
590 591 592
    }


593 594 595 596
};


#endif // FFMMALGORITHMPERIODIC_HPP