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

4 5

#include "../Utils/FGlobal.hpp"
6
#include "../Utils/FGlobalPeriodic.hpp"
7
#include "../Utils/FAssert.hpp"
8
#include "../Utils/FLog.hpp"
9 10
#include "../Utils/FTrace.hpp"
#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>
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(0), OctreeHeight(tree->getHeight()),
BRAMAS Berenger's avatar
BRAMAS Berenger committed
50
          nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3) {
51

52 53
        FAssertLF(tree, "tree cannot be null");
        FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
54

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){
71
        FAssertLF(kernels, "kernels cannot be null");
BRAMAS Berenger's avatar
BRAMAS Berenger committed
72
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
73

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

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

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

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

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

89

90 91 92 93 94 95 96
    /////////////////////////////////////////////////////////////////////////////
    // P2M
    /////////////////////////////////////////////////////////////////////////////

    /** P2M */
    void bottomPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
97
        FLOG( FLog::Controller.write("\tStart Bottom Pass\n").write(FLog::Flush) );
98 99
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
100 101 102 103 104 105 106 107

        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
108
            FLOG(computationCounter.tic());
109
            kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentListSrc());
110
            FLOG(computationCounter.tac());
111 112
        } while(octreeIterator.moveRight());

113 114
        FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
115 116 117 118 119 120 121 122 123
    }

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

    /** M2M */
    void upwardPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
124
        FLOG( FLog::Controller.write("\tStart Upward Pass\n").write(FLog::Flush); );
125 126
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
127 128 129 130 131 132 133 134 135 136

        // 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 ){
137
            FLOG(FTic counterTimeLevel);
138
            const int fackLevel = idxLevel + offsetRealTree;
139 140 141 142
            // for each cells
            do{
                // We need the current cell and the child
                // child is an array (of 8 child) that may be null
143
                FLOG(computationCounter.tic());
144
                kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
145
                FLOG(computationCounter.tac());
146 147 148 149
            } while(octreeIterator.moveRight());

            avoidGotoLeftIterator.moveUp();
            octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
150
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
151 152 153
        }


154 155
        FLOG( FLog::Controller << "\tFinished (@Upward Pass (M2M) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
156 157 158
    }

    /////////////////////////////////////////////////////////////////////////////
159
    // Transfer
160 161 162
    /////////////////////////////////////////////////////////////////////////////

    /** M2L L2L */
163
    void transferPass(){
164 165
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );

166
        FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
167 168
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter);
169 170 171 172

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

173
        const CellClass* neighbors[343];
174 175 176

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

189
            FLOG(computationCounter.tic());
190
            kernels->finishedLevelM2L(fackLevel);
191
            FLOG(computationCounter.tac());
192
            FLOG( FLog::Controller << "\t\t>> Level " << idxLevel << "(" << fackLevel << ") = "  << counterTimeLevel.tacAndElapsed() << "s\n" );
193
        }
194 195
        FLOG( FLog::Controller << "\tFinished (@Downward Pass (M2L) = "  << counterTime.tacAndElapsed() << "s)\n" );
        FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.cumulated() << " s\n" );
196
    }
197

198 199 200 201 202 203 204
    /////////////////////////////////////////////////////////////////////////////
    // Downward
    /////////////////////////////////////////////////////////////////////////////


    void downardPass(){ // second L2L
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
205
        FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
206 207
        FLOG(FTic counterTime);
        FLOG(FTic computationCounter );
208 209 210

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

212 213 214
        const int heightMinusOne = OctreeHeight - 1;
        // for each levels exepted leaf level
        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
215
            FLOG(FTic counterTimeLevel);
216 217
            const int fackLevel = idxLevel + offsetRealTree;

218 219
            // for each cells
            do{
220
                FLOG(computationCounter.tic());
221
                kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), fackLevel);
222
                FLOG(computationCounter.tac());
223 224 225 226
            } while(octreeIterator.moveRight());

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

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

233 234 235 236 237 238 239 240 241 242

    }

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

    /** P2P */
    void directPass(){
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
243
        FLOG( FLog::Controller.write("\tStart Direct Pass\n").write(FLog::Flush); );
244 245 246
        FLOG(FTic counterTime);
        FLOG(FTic computationCounterL2P);
        FLOG(FTic computationCounterP2P);
247 248

        const int heightMinusOne = OctreeHeight - 1;
249
        const FReal boxWidth = tree->getBoxWidth();
250 251 252 253

        typename OctreeClass::Iterator octreeIterator(tree);
        octreeIterator.gotoBottomLeft();
        // There is a maximum of 26 neighbors
berenger-bramas's avatar
berenger-bramas committed
254
        ContainerClass* neighbors[27];
255
        FTreeCoordinate offsets[27];
256
        bool hasPeriodicLeaves;
257 258
        // for each leafs
        do{
259
            FLOG(computationCounterL2P.tic());
260
            kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
261
            FLOG(computationCounterL2P.tac());
262

263
            // need the current particles and neighbors particles
264
            const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
BRAMAS Berenger's avatar
BRAMAS Berenger committed
265
            const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
266 267 268 269 270
            int periodicNeighborsCounter = 0;

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

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
273 274 275 276 277
                    if( neighbors[idxNeig] && !offsets[idxNeig].equals(0,0,0) ){
                        // Put periodic neighbors into other array
                        periodicNeighbors[idxNeig] = neighbors[idxNeig];
                        neighbors[idxNeig] = 0;
                        ++periodicNeighborsCounter;
278 279 280 281 282 283 284 285 286

                        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());
287 288 289 290
                        }
                    }
                }

291
                FLOG(computationCounterP2P.tic());
292 293
                kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                             octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborsCounter);
294
                FLOG(computationCounterP2P.tac());
295 296

                for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
297
                    if( periodicNeighbors[idxNeig] ){
298 299 300 301 302 303 304 305
                        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());
306 307 308 309 310
                        }
                    }
                }
            }

311
            FLOG(computationCounterP2P.tic());
312 313
            kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
                         octreeIterator.getCurrentListSrc(), neighbors, counter - periodicNeighborsCounter);
314
            FLOG(computationCounterP2P.tac());
315 316


317 318 319
        } while(octreeIterator.moveRight());


320 321 322
        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" );
323 324 325 326 327 328 329

    }

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

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 403 404 405 406 407 408 409
    /** 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
410
                    if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1 ){
411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
                        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
441

442
    long long int theoricalRepetition() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
443 444 445 446 447 448 449 450
        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;
451 452
    }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
453

454
    void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
455 456 457 458 459 460 461 462 463 464 465
        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);
        }
466 467 468
    }


469
    FReal extendedBoxWidth() const {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
470 471 472
        // 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));
473 474
    }

475 476 477 478 479 480 481 482
    /** 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
483
        const FReal originalBoxWidth     = tree->getBoxWidth();
484
        const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
485
        const FPoint originalBoxCenter   = tree->getBoxCenter();
486

BRAMAS Berenger's avatar
BRAMAS Berenger committed
487
        const FReal offset = extendedBoxWidth()/FReal(2.0);
488 489 490
        return FPoint( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
                       originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
491 492 493 494 495 496 497 498
    }

    /** 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
      */
499
    int extendedTreeHeight() const {
500 501 502 503
        // The real height
        return OctreeHeight + offsetRealTree;
    }

504 505 506 507 508 509 510 511 512 513 514
    /** 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
      */
515
    void processPeriodicLevels(){
516
        FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
517
        FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
518
        FLOG(FTic counterTime);
519

BRAMAS Berenger's avatar
BRAMAS Berenger committed
520 521 522 523 524 525 526 527
        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);
528
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
529 530 531 532 533
            {
                CellClass* virtualChild[8];
                for(int idxLevel = offsetRealTree-1 ; idxLevel > 1  ; --idxLevel){
                    FMemUtils::setall(virtualChild,&upperCells[idxLevel],8);
                    kernels->M2M( &upperCells[idxLevel-1], virtualChild, idxLevel);
534
                }
535
            }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
536
            CellClass*const downerCells = new CellClass[offsetRealTree];
537

BRAMAS Berenger's avatar
BRAMAS Berenger committed
538 539
            {
                const int idxUpperLevel = 2;
540 541 542 543

                const CellClass* neighbors[343];
                memset(neighbors, 0, sizeof(CellClass*) * 343);
                int counter = 0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
544 545 546 547 548 549
                for(int idxX = -2 ; idxX <= 1 ; ++idxX){
                    for(int idxY = -2 ; idxY <= 1 ; ++idxY){
                        for(int idxZ = -2 ; idxZ <= 1 ; ++idxZ){
                            if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
                                neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
                                ++counter;
550
                            }
551 552 553
                        }
                    }
                }
BRAMAS Berenger's avatar
BRAMAS Berenger committed
554 555 556
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
557

BRAMAS Berenger's avatar
BRAMAS Berenger committed
558 559 560 561 562 563 564 565 566 567
            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;
568 569 570
                            }
                        }
                    }
571 572
                }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
573 574 575
                // compute M2L
                kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
            }
576

BRAMAS Berenger's avatar
BRAMAS Berenger committed
577
            {
578 579
                CellClass* virtualChild[8];
                memset(virtualChild, 0, sizeof(CellClass*) * 8);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
580 581 582 583
                for(int idxLevel = 2 ; idxLevel <= offsetRealTree-1  ; ++idxLevel){
                    virtualChild[0] = &downerCells[idxLevel];
                    kernels->L2L( &downerCells[idxLevel-1], virtualChild, idxLevel);
                }
584
            }
585

BRAMAS Berenger's avatar
BRAMAS Berenger committed
586 587 588 589 590
            // L2L from 0 to level 1
            {
                typename OctreeClass::Iterator octreeIterator(tree);
                octreeIterator.gotoLeft();
                kernels->L2L( &downerCells[offsetRealTree-1], octreeIterator.getCurrentBox(), offsetRealTree);
591 592
            }

BRAMAS Berenger's avatar
BRAMAS Berenger committed
593 594
            delete[] upperCells;
            delete[] downerCells;
595
        }
596

597
        FLOG( FLog::Controller << "\tFinished (@Periodic = "  << counterTime.tacAndElapsed() << "s)\n" );
598 599 600
    }


601 602 603 604
};


#endif // FFMMALGORITHMPERIODIC_HPP