FGroupTree.hpp 65.6 KB
Newer Older
1 2

// Keep in private GIT
3 4
#ifndef FGROUPTREE_HPP
#define FGROUPTREE_HPP
BRAMAS Berenger's avatar
BRAMAS Berenger committed
5
#include <vector>
6
#include <functional>
7

8 9 10 11
#include "../../Utils/FAssert.hpp"
#include "../../Utils/FPoint.hpp"
#include "../../Utils/FQuickSort.hpp"
#include "../../Containers/FTreeCoordinate.hpp"
12
#include "../../Containers/FCoordinateComputer.hpp"
13 14
#include "FGroupOfCells.hpp"
#include "FGroupOfParticles.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
15
#include "FGroupAttachedLeaf.hpp"
16
#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"
17

18

19

20
template <class FReal, class SymbolCellClass, class PoleCellClass, class LocalCellClass,
21
          class GroupAttachedLeafClass, unsigned NbSymbAttributes, unsigned NbAttributesPerParticle, class AttributeClass = FReal>
22
class FGroupTree {
23
public:
24
    typedef GroupAttachedLeafClass BasicAttachedClass;
25
    typedef FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle,AttributeClass> ParticleGroupClass;
26
    typedef FGroupOfCells<SymbolCellClass, PoleCellClass, LocalCellClass> CellGroupClass;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
27

28
protected:
29
    //< height of the tree (1 => only the root)
30
    const int _treeHeight;
31
    //< max number of cells in a block
32
    const int _nbElementsPerBlock;
33
    //< all the blocks of the tree
34
    std::vector<CellGroupClass*>* _cellBlocksPerLevel;
35
    //< all the blocks of leaves
36
    std::vector<ParticleGroupClass*> _particleBlocks;
37

BRAMAS Berenger's avatar
BRAMAS Berenger committed
38
    //< the space system center
39
    const FPoint<FReal> boxCenter;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
40
    //< the space system corner (used to compute morton index)
41
    const FPoint<FReal> boxCorner;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
42 43 44 45 46
    //< the space system width
    const FReal boxWidth;
    //< the width of a box at width level
    const FReal boxWidthAtLeafLevel;

47
public:
BRAMAS Berenger's avatar
BRAMAS Berenger committed
48 49 50 51 52
    typedef typename std::vector<CellGroupClass*>::iterator CellGroupIterator;
    typedef typename std::vector<CellGroupClass*>::const_iterator CellGroupConstIterator;
    typedef typename std::vector<ParticleGroupClass*>::iterator ParticleGroupIterator;
    typedef typename std::vector<ParticleGroupClass*>::const_iterator ParticleGroupConstIterator;

53
    /** This constructor create a blocked octree from a usual octree
54 55 56 57
     * The cell are allocated as in the usual octree (no copy constructor are called!)
     * Once allocated each cell receive its morton index and tree coordinate.
     * No blocks are allocated at level 0.
     */
58
    template<class OctreeClass>
59
    FGroupTree(const int in_treeHeight, const int in_nbElementsPerBlock, OctreeClass*const inOctreeSrc)
60
        : _treeHeight(in_treeHeight), _nbElementsPerBlock(in_nbElementsPerBlock), _cellBlocksPerLevel(nullptr),
BRAMAS Berenger's avatar
BRAMAS Berenger committed
61
          boxCenter(inOctreeSrc->getBoxCenter()), boxCorner(inOctreeSrc->getBoxCenter(),-(inOctreeSrc->getBoxWidth()/2)),
62
          boxWidth(inOctreeSrc->getBoxWidth()), boxWidthAtLeafLevel(inOctreeSrc->getBoxWidth()/FReal(1<<(in_treeHeight-1))){
63
        _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];
64 65 66 67 68 69

        // Iterate on the tree and build
        typename OctreeClass::Iterator octreeIterator(inOctreeSrc);
        octreeIterator.gotoBottomLeft();

        { // First leaf level, we create leaves and cells groups
70
            const int idxLevel = _treeHeight-1;
71 72 73 74
            typename OctreeClass::Iterator avoidGotoLeft = octreeIterator;
            // For each cell at this level
            do {
                typename OctreeClass::Iterator blockIteratorInOctree = octreeIterator;
75
                // Move the iterator per _nbElementsPerBlock (or until it cannot move right)
76
                int sizeOfBlock = 1;
77
                FSize nbParticlesInGroup = octreeIterator.getCurrentLeaf()->getSrc()->getNbParticles();
78
                while(sizeOfBlock < _nbElementsPerBlock && octreeIterator.moveRight()){
79 80 81 82 83
                    sizeOfBlock += 1;
                    nbParticlesInGroup += octreeIterator.getCurrentLeaf()->getSrc()->getNbParticles();
                }

                // Create a block with the apropriate parameters
84
                CellGroupClass*const newBlock = new CellGroupClass(blockIteratorInOctree.getCurrentGlobalIndex(),
85 86
                                                                 octreeIterator.getCurrentGlobalIndex()+1,
                                                                 sizeOfBlock);
87
                FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>(blockIteratorInOctree.getCurrentGlobalIndex(),
88 89
                                octreeIterator.getCurrentGlobalIndex()+1,
                                 sizeOfBlock, nbParticlesInGroup);
90 91 92

                // Initialize each cell of the block
                int cellIdInBlock = 0;
93
                size_t nbParticlesOffsetBeforeLeaf = 0;
94
                while(cellIdInBlock != sizeOfBlock){
95 96
                    const MortonIndex newNodeIndex = blockIteratorInOctree.getCurrentCell()->getMortonIndex();
                    const FTreeCoordinate newNodeCoordinate = blockIteratorInOctree.getCurrentCell()->getCoordinate();
97
                    // Add cell
98
                    newBlock->newCell(newNodeIndex, cellIdInBlock);
99

100 101 102 103
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(newNodeIndex);
                    symbolic.setCoordinate(newNodeCoordinate);
                    symbolic.setLevel(idxLevel);
104 105

                    // Add leaf
106
                    nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(newNodeIndex, cellIdInBlock,
107
                                              blockIteratorInOctree.getCurrentLeaf()->getSrc()->getNbParticles(),
108
                                              nbParticlesOffsetBeforeLeaf);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
109

110
                    BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(cellIdInBlock);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
111
                    attachedLeaf.copyFromContainer(blockIteratorInOctree.getCurrentLeaf()->getSrc(), 0);
112 113 114 115 116 117

                    cellIdInBlock += 1;
                    blockIteratorInOctree.moveRight();
                }

                // Keep the block
118 119
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
                _particleBlocks.push_back(newParticleBlock);
120 121 122 123 124 125 126 127 128

                // If we can move right then add another block
            } while(octreeIterator.moveRight());

            avoidGotoLeft.moveUp();
            octreeIterator = avoidGotoLeft;
        }

        // For each level from heigth - 2 to 1
129
        for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
130 131 132 133
            typename OctreeClass::Iterator avoidGotoLeft = octreeIterator;
            // For each cell at this level
            do {
                typename OctreeClass::Iterator blockIteratorInOctree = octreeIterator;
134
                // Move the iterator per _nbElementsPerBlock (or until it cannot move right)
135
                int sizeOfBlock = 1;
136
                while(sizeOfBlock < _nbElementsPerBlock && octreeIterator.moveRight()){
137 138 139 140
                    sizeOfBlock += 1;
                }

                // Create a block with the apropriate parameters
141
                CellGroupClass*const newBlock = new CellGroupClass(blockIteratorInOctree.getCurrentGlobalIndex(),
142 143 144 145 146
                                                                 octreeIterator.getCurrentGlobalIndex()+1,
                                                                 sizeOfBlock);
                // Initialize each cell of the block
                int cellIdInBlock = 0;
                while(cellIdInBlock != sizeOfBlock){
147 148 149
                    const MortonIndex newNodeIndex = blockIteratorInOctree.getCurrentCell()->getMortonIndex();
                    const FTreeCoordinate newNodeCoordinate = blockIteratorInOctree.getCurrentCell()->getCoordinate();
                    newBlock->newCell(newNodeIndex, cellIdInBlock);
150

151 152 153 154
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(newNodeIndex);
                    symbolic.setCoordinate(newNodeCoordinate);
                    symbolic.setLevel(idxLevel);
155 156 157 158 159 160

                    cellIdInBlock += 1;
                    blockIteratorInOctree.moveRight();
                }

                // Keep the block
161
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
162 163 164 165 166 167 168 169 170

                // If we can move right then add another block
            } while(octreeIterator.moveRight());

            avoidGotoLeft.moveUp();
            octreeIterator = avoidGotoLeft;
        }
    }

171 172 173 174 175 176 177 178 179
    /**
     * This constructor create a group tree from a particle container index.
     * The morton index are computed and the particles are sorted in a first stage.
     * Then the leaf level is done.
     * Finally the other leve are proceed one after the other.
     * It should be easy to make it parallel using for and tasks.
     * If no limite give inLeftLimite = -1
     */
    template<class ParticleContainer>
180 181
    FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
               const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
182
               const bool particlesAreSorted = false, MortonIndex inLeftLimite = -1):
183
            _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
184
            boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
185 186
            boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
    {
187

188
        _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];
189

190
        MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
191 192 193 194
        // First we work at leaf level
        {
            // Build morton index for particles
            struct ParticleSortingStruct{
195
                FSize originalIndex;
196 197 198
                MortonIndex mindex;
            };
            // Convert position to morton index
199
            const FSize nbParticles = inParticlesContainer->getNbParticles();
200 201 202 203 204 205
            ParticleSortingStruct* particlesToSort = new ParticleSortingStruct[nbParticles];
            {
                const FReal* xpos = inParticlesContainer->getPositions()[0];
                const FReal* ypos = inParticlesContainer->getPositions()[1];
                const FReal* zpos = inParticlesContainer->getPositions()[2];

206
                for(FSize idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
207
                    const FTreeCoordinate host = FCoordinateComputer::GetCoordinateFromPositionAndCorner<FReal>(this->boxCorner, this->boxWidth,
208
                                                                                                       _treeHeight,
209
                                                                                                       FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
210
                    const MortonIndex particleIndex = host.getMortonIndex();
211 212 213 214 215 216 217
                    particlesToSort[idxPart].mindex = particleIndex;
                    particlesToSort[idxPart].originalIndex = idxPart;
                }
            }

            // Sort if needed
            if(particlesAreSorted == false){
218
                FQuickSort<ParticleSortingStruct, FSize>::QsOmp(particlesToSort, nbParticles, [](const ParticleSortingStruct& v1, const ParticleSortingStruct& v2){
219 220 221 222 223 224
                    return v1.mindex <= v2.mindex;
                });
            }

            FAssertLF(nbParticles == 0 || inLeftLimite < particlesToSort[0].mindex);
            // Convert to block
225 226
            const int idxLevel = (_treeHeight - 1);
            FSize* nbParticlesPerLeaf = new FSize[_nbElementsPerBlock];
227
            FSize firstParticle = 0;
228 229 230
            // We need to proceed each group in sub level
            while(firstParticle != nbParticles){
                int sizeOfBlock = 0;
231
                FSize lastParticle = firstParticle;
232
                // Count until end of sub group is reached or we have enough cells
233
                while(sizeOfBlock < _nbElementsPerBlock && lastParticle < nbParticles){
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
                    if(sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != particlesToSort[lastParticle].mindex){
                        currentBlockIndexes[sizeOfBlock] = particlesToSort[lastParticle].mindex;
                        nbParticlesPerLeaf[sizeOfBlock]  = 1;
                        sizeOfBlock += 1;
                    }
                    else{
                        nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    }
                    lastParticle += 1;
                }
                while(lastParticle < nbParticles && currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex){
                    nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    lastParticle += 1;
                }

                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
251 252
                               currentBlockIndexes[sizeOfBlock-1]+1,
                               sizeOfBlock);
253
                FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>(currentBlockIndexes[0],
254 255 256 257 258
                        currentBlockIndexes[sizeOfBlock-1]+1,
                        sizeOfBlock, lastParticle-firstParticle);

                // Init cells
                size_t nbParticlesOffsetBeforeLeaf = 0;
259
                FSize offsetParticles = firstParticle;
260 261 262
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

263 264
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
265
                    FTreeCoordinate coord;
266
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
267 268
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
269 270 271 272 273

                    // Add leaf
                    nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(currentBlockIndexes[cellIdInBlock], cellIdInBlock,
                                              nbParticlesPerLeaf[cellIdInBlock], nbParticlesOffsetBeforeLeaf);

274
                    BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(cellIdInBlock);
275
                    // Copy each particle from the original position
276
                    for(FSize idxPart = 0 ; idxPart < nbParticlesPerLeaf[cellIdInBlock] ; ++idxPart){
277 278 279 280 281 282
                        attachedLeaf.setParticle(idxPart, particlesToSort[idxPart + offsetParticles].originalIndex, inParticlesContainer);
                    }
                    offsetParticles += nbParticlesPerLeaf[cellIdInBlock];
                }

                // Keep the block
283 284
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
                _particleBlocks.push_back(newParticleBlock);
285 286 287 288 289 290 291 292 293

                sizeOfBlock = 0;
                firstParticle = lastParticle;
            }
            delete[] nbParticlesPerLeaf;
            delete[] particlesToSort;
        }

        // For each level from heigth - 2 to 1
294
        for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
295 296
            inLeftLimite = (inLeftLimite == -1 ? inLeftLimite : (inLeftLimite>>3));

297 298
            CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
            const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316

            // Skip blocks that do not respect limit
            while(iterChildCells != iterChildEndCells
                  && ((*iterChildCells)->getEndingIndex()>>3) <= inLeftLimite){
                ++iterChildCells;
            }
            // If lower level is empty or all blocks skiped stop here
            if(iterChildCells == iterChildEndCells){
                break;
            }

            MortonIndex currentCellIndex = (*iterChildCells)->getStartingIndex();
            if((currentCellIndex>>3) <= inLeftLimite) currentCellIndex = ((inLeftLimite+1)<<3);
            int sizeOfBlock = 0;

            // We need to proceed each group in sub level
            while(iterChildCells != iterChildEndCells){
                // Count until end of sub group is reached or we have enough cells
317
                while(sizeOfBlock < _nbElementsPerBlock && iterChildCells != iterChildEndCells ){
318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
                    if((sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != (currentCellIndex>>3))
                            && (*iterChildCells)->exists(currentCellIndex)){
                        currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                        sizeOfBlock += 1;
                        currentCellIndex = (((currentCellIndex>>3)+1)<<3);
                    }
                    else{
                        currentCellIndex += 1;
                    }
                    // If we are at the end of the sub group, move to next
                    while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
                        ++iterChildCells;
                        // Update morton index
                        if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                            currentCellIndex = (*iterChildCells)->getStartingIndex();
                        }
                    }
                }

                // If group is full
338
                if(sizeOfBlock == _nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
339 340
                    // Create a group
                    CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
341 342
                                    currentBlockIndexes[sizeOfBlock-1]+1,
                                    sizeOfBlock);
343 344 345 346
                    // Init cells
                    for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                        newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

347 348
                        SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                        symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
349
                        FTreeCoordinate coord;
350
                        coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
351 352
                        symbolic.setCoordinate(coord);
                        symbolic.setLevel(idxLevel);
353 354 355
                    }

                    // Keep the block
356
                    _cellBlocksPerLevel[idxLevel].push_back(newBlock);
357 358 359 360 361 362 363 364

                    sizeOfBlock = 0;
                }
            }
        }
        delete[] currentBlockIndexes;
    }

365 366 367 368 369 370 371 372 373 374 375 376 377 378
    /**
     * This constructor create a group tree from a particle container index.
     * The morton index are computed and the particles are sorted in a first stage.
     * Then the leaf level is done.
     * Finally the other leve are proceed one after the other.
     * It should be easy to make it parallel using for and tasks.
     * If no limite give inLeftLimite = -1
     * The cover ration is the minimum pourcentage of cell that should
     * exist in a group (0 means no limite, 1 means the block must be dense)
     * oneParent should be turned on if it is better to have one block parent
     * per sublock (in case of have the cost of FMM that increase with the level
     * this could be an asset).
     */
    template<class ParticleContainer>
379 380
    FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
               const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
381 382
               const bool particlesAreSorted, const bool oneParent,
               const FReal inCoverRatio = 0.0, MortonIndex inLeftLimite = -1):
383
            _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
384
            boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
385 386
            boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
    {
387

388 389 390
        FAssertLF(inCoverRatio == 0.0 || oneParent == true, "If a ratio is choosen oneParent should be turned on");
        const bool userCoverRatio = (inCoverRatio != 0.0);

391
        _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];
392

393
        MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
394 395 396 397
        // First we work at leaf level
        {
            // Build morton index for particles
            struct ParticleSortingStruct{
398
                FSize originalIndex;
399 400 401
                MortonIndex mindex;
            };
            // Convert position to morton index
402
            const FSize nbParticles = inParticlesContainer->getNbParticles();
403 404 405 406 407 408
            ParticleSortingStruct* particlesToSort = new ParticleSortingStruct[nbParticles];
            {
                const FReal* xpos = inParticlesContainer->getPositions()[0];
                const FReal* ypos = inParticlesContainer->getPositions()[1];
                const FReal* zpos = inParticlesContainer->getPositions()[2];

409
                for(FSize idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
410
                    const FTreeCoordinate host = FCoordinateComputer::GetCoordinateFromPositionAndCorner<FReal>(this->boxCorner, this->boxWidth,
411
                                                                                                       _treeHeight,
412
                                                                                                       FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
413
                    const MortonIndex particleIndex = host.getMortonIndex();
414 415 416 417 418 419 420
                    particlesToSort[idxPart].mindex = particleIndex;
                    particlesToSort[idxPart].originalIndex = idxPart;
                }
            }

            // Sort if needed
            if(particlesAreSorted == false){
421
                FQuickSort<ParticleSortingStruct, FSize>::QsOmp(particlesToSort, nbParticles, [](const ParticleSortingStruct& v1, const ParticleSortingStruct& v2){
422 423 424 425 426 427 428
                    return v1.mindex <= v2.mindex;
                });
            }

            FAssertLF(nbParticles == 0 || inLeftLimite < particlesToSort[0].mindex);

            // Convert to block
429 430
            const int idxLevel = (_treeHeight - 1);
            int* nbParticlesPerLeaf = new int[_nbElementsPerBlock];
431 432 433 434 435 436
            int firstParticle = 0;
            // We need to proceed each group in sub level
            while(firstParticle != nbParticles){
                int sizeOfBlock = 0;
                int lastParticle = firstParticle;
                // Count until end of sub group is reached or we have enough cells
437
                while(sizeOfBlock < _nbElementsPerBlock && lastParticle < nbParticles
438 439 440 441
                      && (userCoverRatio == false
                          || sizeOfBlock == 0
                          || currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex
                          || (FReal(sizeOfBlock+1)/FReal(particlesToSort[lastParticle].mindex-particlesToSort[firstParticle].mindex)) >= inCoverRatio)){
442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470
                    if(sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != particlesToSort[lastParticle].mindex){
                        currentBlockIndexes[sizeOfBlock] = particlesToSort[lastParticle].mindex;
                        nbParticlesPerLeaf[sizeOfBlock]  = 1;
                        sizeOfBlock += 1;
                    }
                    else{
                        nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    }
                    lastParticle += 1;
                }
                while(lastParticle < nbParticles && currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex){
                    nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    lastParticle += 1;
                }

                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
                                                                 currentBlockIndexes[sizeOfBlock-1]+1,
                                                                 sizeOfBlock);
                FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>(currentBlockIndexes[0],
                        currentBlockIndexes[sizeOfBlock-1]+1,
                        sizeOfBlock, lastParticle-firstParticle);

                // Init cells
                size_t nbParticlesOffsetBeforeLeaf = 0;
                int offsetParticles = firstParticle;
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

471 472
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
473
                    FTreeCoordinate coord;
474
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
475 476
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
477 478 479 480 481

                    // Add leaf
                    nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(currentBlockIndexes[cellIdInBlock], cellIdInBlock,
                                              nbParticlesPerLeaf[cellIdInBlock], nbParticlesOffsetBeforeLeaf);

482
                    BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(cellIdInBlock);
483
                    // Copy each particle from the original position
484
                    for(FSize idxPart = 0 ; idxPart < nbParticlesPerLeaf[cellIdInBlock] ; ++idxPart){
485 486 487 488 489 490
                        attachedLeaf.setParticle(idxPart, particlesToSort[idxPart + offsetParticles].originalIndex, inParticlesContainer);
                    }
                    offsetParticles += nbParticlesPerLeaf[cellIdInBlock];
                }

                // Keep the block
491 492
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
                _particleBlocks.push_back(newParticleBlock);
493 494 495 496 497 498 499 500 501 502

                sizeOfBlock = 0;
                firstParticle = lastParticle;
            }
            delete[] nbParticlesPerLeaf;
            delete[] particlesToSort;
        }


        // For each level from heigth - 2 to 1
503
        for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
504 505
            inLeftLimite = (inLeftLimite == -1 ? inLeftLimite : (inLeftLimite>>3));

506 507
            CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
            const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526

            // Skip blocks that do not respect limit
            while(iterChildCells != iterChildEndCells
                  && ((*iterChildCells)->getEndingIndex()>>3) <= inLeftLimite){
                ++iterChildCells;
            }
            // If lower level is empty or all blocks skiped stop here
            if(iterChildCells == iterChildEndCells){
                break;
            }

            MortonIndex currentCellIndex = (*iterChildCells)->getStartingIndex();
            if((currentCellIndex>>3) <= inLeftLimite) currentCellIndex = ((inLeftLimite+1)<<3);
            int sizeOfBlock = 0;

            if(oneParent == false){
                // We need to proceed each group in sub level
                while(iterChildCells != iterChildEndCells){
                    // Count until end of sub group is reached or we have enough cells
527
                    while(sizeOfBlock < _nbElementsPerBlock && iterChildCells != iterChildEndCells ){
528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
                        if((sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != (currentCellIndex>>3))
                                && (*iterChildCells)->exists(currentCellIndex)){
                            currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                            sizeOfBlock += 1;
                            currentCellIndex = (((currentCellIndex>>3)+1)<<3);
                        }
                        else{
                            currentCellIndex += 1;
                        }
                        // If we are at the end of the sub group, move to next
                        while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
                            ++iterChildCells;
                            // Update morton index
                            if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                                currentCellIndex = (*iterChildCells)->getStartingIndex();
                            }
                        }
                    }

                    // If group is full
548
                    if(sizeOfBlock == _nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
549 550 551 552 553 554 555 556
                        // Create a group
                        CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
                                                                         currentBlockIndexes[sizeOfBlock-1]+1,
                                                                         sizeOfBlock);
                        // Init cells
                        for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                            newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

557 558
                            SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                            symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
559
                            FTreeCoordinate coord;
560
                            coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
561 562
                            symbolic.setCoordinate(coord);
                            symbolic.setLevel(idxLevel);
563 564 565
                        }

                        // Keep the block
566
                        _cellBlocksPerLevel[idxLevel].push_back(newBlock);
567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593

                        sizeOfBlock = 0;
                    }
                }
            }
            else{
                // We need to proceed each group in sub level
                while(iterChildCells != iterChildEndCells){
                    // We want one parent group per child group so we will stop the parent group
                    // when we arrive to the same parent as lastChildIndex (which is lastChildIndex>>3)
                    const MortonIndex lastChildIndex = ((*iterChildCells)->getEndingIndex()-1);
                    // Count until end of sub group is reached or we passe the requested parent
                    while( iterChildCells != iterChildEndCells
                           && (currentCellIndex>>3) <= (lastChildIndex>>3) ){
                        // Proceed until the requested parent
                        while(currentCellIndex != (*iterChildCells)->getEndingIndex()
                              && (currentCellIndex>>3) <= (lastChildIndex>>3) ){
                            if((*iterChildCells)->exists(currentCellIndex)){
                                currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                                sizeOfBlock += 1;
                                currentCellIndex = (((currentCellIndex>>3)+1)<<3);
                            }
                            else{
                                currentCellIndex += 1;
                            }
                        }
                        // If we are at the end of the sub group, move to next (otherwise we have consume a part of it)
594
                        while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
                            ++iterChildCells;
                            // Update morton index
                            if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                                currentCellIndex = (*iterChildCells)->getStartingIndex();
                            }
                        }
                    }

                    // If group is full
                    if(sizeOfBlock){
                        // Create a group
                        CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
                                                                         currentBlockIndexes[sizeOfBlock-1]+1,
                                                                         sizeOfBlock);
                        // Init cells
                        for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                            newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

613 614
                            SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                            symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
615
                            FTreeCoordinate coord;
616
                            coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
617 618
                            symbolic.setCoordinate(coord);
                            symbolic.setLevel(idxLevel);
619 620 621
                        }

                        // Keep the block
622
                        _cellBlocksPerLevel[idxLevel].push_back(newBlock);
623 624 625 626 627 628 629 630 631

                        sizeOfBlock = 0;
                    }
                }
            }
        }
        delete[] currentBlockIndexes;
    }

632
    template<class ParticleContainer>
633 634
    FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
               const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
635
			   std::vector<std::vector<int>> & blockSizeAtEachLevel,
636
               const bool particlesAreSorted = false):
637
            _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
638
            boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
639 640
            boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
    {
641

642
        _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];
643

644
        MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
        // First we work at leaf level
        {
            // Build morton index for particles
            struct ParticleSortingStruct{
                FSize originalIndex;
                MortonIndex mindex;
            };
            // Convert position to morton index
            const FSize nbParticles = inParticlesContainer->getNbParticles();
            ParticleSortingStruct* particlesToSort = new ParticleSortingStruct[nbParticles];
            {
                const FReal* xpos = inParticlesContainer->getPositions()[0];
                const FReal* ypos = inParticlesContainer->getPositions()[1];
                const FReal* zpos = inParticlesContainer->getPositions()[2];

                for(FSize idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
                    const FTreeCoordinate host = FCoordinateComputer::GetCoordinateFromPositionAndCorner<FReal>(this->boxCorner, this->boxWidth,
662
                                                                                                       _treeHeight,
663
                                                                                                       FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
664
                    const MortonIndex particleIndex = host.getMortonIndex();
665 666 667 668 669 670 671 672 673 674 675 676 677 678
                    particlesToSort[idxPart].mindex = particleIndex;
                    particlesToSort[idxPart].originalIndex = idxPart;
                }
            }

            // Sort if needed
            if(particlesAreSorted == false){
                FQuickSort<ParticleSortingStruct, FSize>::QsOmp(particlesToSort, nbParticles, [](const ParticleSortingStruct& v1, const ParticleSortingStruct& v2){
                    return v1.mindex <= v2.mindex;
                });
            }


            // Convert to block
679
            const int idxLevel = (_treeHeight - 1);
680
			int idxBlock = 0;
681
            FSize* nbParticlesPerLeaf = new FSize[_nbElementsPerBlock];
682 683 684 685 686 687
            FSize firstParticle = 0;
            // We need to proceed each group in sub level
            while(firstParticle != nbParticles){
                int sizeOfBlock = 0;
                FSize lastParticle = firstParticle;
                // Count until end of sub group is reached or we have enough cells
688
                while(sizeOfBlock < blockSizeAtEachLevel[_treeHeight-1][idxBlock] && lastParticle < nbParticles){
689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710
                    if(sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != particlesToSort[lastParticle].mindex){
                        currentBlockIndexes[sizeOfBlock] = particlesToSort[lastParticle].mindex;
                        nbParticlesPerLeaf[sizeOfBlock]  = 1;
                        sizeOfBlock += 1;
                    }
                    else{
                        nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    }
                    lastParticle += 1;
                }
                while(lastParticle < nbParticles && currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex){
                    nbParticlesPerLeaf[sizeOfBlock-1] += 1;
                    lastParticle += 1;
                }

                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
                                                                 currentBlockIndexes[sizeOfBlock-1]+1,
                                                                 sizeOfBlock);
                FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>(currentBlockIndexes[0],
                        currentBlockIndexes[sizeOfBlock-1]+1,
                        sizeOfBlock, lastParticle-firstParticle);
711
                        #include <iostream>
712 713 714
				using namespace std;
				if(currentBlockIndexes[sizeOfBlock-1]+1 == 511)
					cout << "Suricate" << endl;
715 716 717 718 719 720
                // Init cells
                size_t nbParticlesOffsetBeforeLeaf = 0;
                FSize offsetParticles = firstParticle;
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

721 722
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
723
                    FTreeCoordinate coord;
724 725 726
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
727 728 729 730 731 732 733 734 735 736 737 738 739 740

                    // Add leaf
                    nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(currentBlockIndexes[cellIdInBlock], cellIdInBlock,
                                              nbParticlesPerLeaf[cellIdInBlock], nbParticlesOffsetBeforeLeaf);

                    BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(cellIdInBlock);
                    // Copy each particle from the original position
                    for(FSize idxPart = 0 ; idxPart < nbParticlesPerLeaf[cellIdInBlock] ; ++idxPart){
                        attachedLeaf.setParticle(idxPart, particlesToSort[idxPart + offsetParticles].originalIndex, inParticlesContainer);
                    }
                    offsetParticles += nbParticlesPerLeaf[cellIdInBlock];
                }

                // Keep the block
741 742
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
                _particleBlocks.push_back(newParticleBlock);
743 744 745

                sizeOfBlock = 0;
                firstParticle = lastParticle;
746
				++idxBlock;
747 748 749 750 751 752 753
            }
            delete[] nbParticlesPerLeaf;
            delete[] particlesToSort;
        }


        // For each level from heigth - 2 to 1
754
        for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
755

756 757
            CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
            const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
758 759 760 761 762 763 764 765

            // If lower level is empty or all blocks skiped stop here
            if(iterChildCells == iterChildEndCells){
                break;
            }

            MortonIndex currentCellIndex = (*iterChildCells)->getStartingIndex();
            int sizeOfBlock = 0;
766
			int idxBlock = 0;
767 768 769
            // We need to proceed each group in sub level
            while(iterChildCells != iterChildEndCells){
                // Count until end of sub group is reached or we have enough cells
770
                while(sizeOfBlock < blockSizeAtEachLevel[idxLevel][idxBlock] && iterChildCells != iterChildEndCells ){
771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790
                    if((sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != (currentCellIndex>>3))
                            && (*iterChildCells)->exists(currentCellIndex)){
                        currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                        sizeOfBlock += 1;
                        currentCellIndex = (((currentCellIndex>>3)+1)<<3);
                    }
                    else{
                        currentCellIndex += 1;
                    }
                    // If we are at the end of the sub group, move to next
                    while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
                        ++iterChildCells;
                        // Update morton index
                        if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                            currentCellIndex = (*iterChildCells)->getStartingIndex();
                        }
                    }
                }

                // If group is full
791
                if(sizeOfBlock == blockSizeAtEachLevel[idxLevel][idxBlock] || (sizeOfBlock && iterChildCells == iterChildEndCells)){ //NOTE la seconde partie va sûrement sauter, car la taille est pré-calculée
792 793 794 795 796 797 798 799
                    // Create a group
                    CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
                                                                     currentBlockIndexes[sizeOfBlock-1]+1,
                                                                     sizeOfBlock);
                    // Init cells
                    for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                        newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

800 801
                        SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                        symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
802
                        FTreeCoordinate coord;
803 804 805
                        coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
                        symbolic.setCoordinate(coord);
                        symbolic.setLevel(idxLevel);
806 807 808
                    }

                    // Keep the block
809
                    _cellBlocksPerLevel[idxLevel].push_back(newBlock);
810 811

                    sizeOfBlock = 0;
812
					++idxBlock;
813 814 815 816 817
                }
            }
        }
        delete[] currentBlockIndexes;
    }
818

819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844

    /**
     * Minimal Constructor of GroupTree
     * @author benjamin.dufoyer@inria.fr
     * @param in__treeHeight size of the tree
     * @param in__nbElementsPerBlock block size
     * @param in_boxCenter  box center
     * @param in_boxCorner  box cornet
     * @param in_boxWidth   bow witdh
     * @param in_boxWidthAtLeafLevel  box width at leaf level
     */
    FGroupTree(
         int in__treeHeight,
         int in__nbElementsPerBlock,
         FPoint<FReal> in_boxCenter,
         FPoint<FReal> in_boxCorner,
         FReal in_boxWidth,
         FReal in_boxWidthAtLeafLevel
    ):
        _treeHeight(in__treeHeight),
        _nbElementsPerBlock(in__nbElementsPerBlock),
        boxCenter(in_boxCenter),
        boxCorner(in_boxCorner),
        boxWidth(in_boxWidth),
        boxWidthAtLeafLevel(in_boxWidthAtLeafLevel)
    {
845
        this->_cellBlocksPerLevel = new std::vector<CellGroupClass*>[this->_treeHeight];
846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894
    }


    /**
     * get_block_tree_instance return a new instance of FGroupTree from
     * a blocked linear tree
     * @author benjamin.dufoyer@inria.fr
     * @param  blocked_linear_tree blocked linear tree
     * @return new FGroupTree
     */
    template<
    class GroupCellSymbClass,
    class GroupCellUpClass,
    class GroupCellDownClass,
    class GroupContainerClass
    >
    static FGroupTree get_block_tree_instance(
        int in_tree_height,
        int in_block_size,
        FPoint<FReal> in_box_center,
        FReal in_box_width
    ){
        // Compute every information to initialise the group tree
        FPoint<FReal> box_corner = FPoint<FReal>(in_box_center, -in_box_width/2);
        FReal box_width_at_leaf_level = in_box_width/FReal( 1<< (in_tree_height-1));

        // Return a new instance of a empty group tree
        return FGroupTree<FReal,GroupCellSymbClass,GroupCellUpClass, GroupCellDownClass, GroupContainerClass, NbSymbAttributes, NbAttributesPerParticle, FReal>(
                 in_tree_height
                ,in_block_size
                ,in_box_center
                ,box_corner
                ,in_box_width
                ,box_width_at_leaf_level);
    }


    /////////////////////////////////////////////////////////
    // Function to init group tree
    /////////////////////////////////////////////////////////

    /**
     * create_tree this function fill the tree from blocked_linear_tree
     * She build the group tree from the bottom
     * @author benjamin.dufoyer@inria.fr
     * @param  in_lin_tree  Blocked linear tree
     * @param  particles    vector where particle are stock,
     *                      they will be sort BEFORE calling this function
     */
895
    template<class Group_Linear_tree,
896 897
              class Particle_Container
            >
898
    void create_tree(Group_Linear_tree* in_lin_tree,
899 900 901 902 903 904 905 906 907
                     Particle_Container* particles){
        // Creation of the leaf level and groups of particle
        create_leaf_level(in_lin_tree,particles);
        // Creation of every level of the tree
        for(int level = this->_treeHeight-2 ; level >= 0 ; --level){
            create_block_nodes_at_level(level);
        }
    }

908 909
    /** This function dealloc the tree by deleting each block */
    ~FGroupTree(){
910
        for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
911
            std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
912
            for (CellGroupClass* block: levelBlocks){
913 914 915
                delete block;
            }
        }
916
        delete[] _cellBlocksPerLevel;
917

918
        for (ParticleGroupClass* block: _particleBlocks){
919 920 921 922
            delete block;
        }
    }

923
    ////////////////////////////////////////////////////////
924 925 926 927 928 929 930 931 932
    // Lambda function to apply to all member
    /////////////////////////////////////////////////////////

    /**
   * @brief forEachLeaf iterate on the leaf and apply the function
   * @param function
   */
    template<class ParticlesAttachedClass>
    void forEachLeaf(std::function<void(ParticlesAttachedClass*)> function){
933
        for (ParticleGroupClass* block: _particleBlocks){
934 935 936 937 938 939 940 941
            block->forEachLeaf(function);
        }
    }

    /**
   * @brief forEachLeaf iterate on the cell and apply the function
   * @param function
   */
942
    void forEachCell(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*)> function){
943
        for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
944
            std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
945
            for (CellGroupClass* block: levelBlocks){
946 947 948 949 950 951 952 953 954
                block->forEachCell(function);
            }
        }
    }

    /**
   * @brief forEachLeaf iterate on the cell and apply the function
   * @param function
   */
955
    void forEachCellWithLevel(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,const int)> function){
956
        for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
957
            std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
958
            for (CellGroupClass* block: levelBlocks){
959 960 961 962 963 964 965 966 967 968
                block->forEachCell(function, idxLevel);
            }
        }
    }

    /**
   * @brief forEachLeaf iterate on the cell and apply the function
   * @param function
   */
    template<class ParticlesAttachedClass>
969
    void forEachCellLeaf(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,ParticlesAttachedClass*)> function){
970 971
        CellGroupIterator iterCells = _cellBlocksPerLevel[_treeHeight-1].begin();
        const CellGroupIterator iterEndCells = _cellBlocksPerLevel[_treeHeight-1].end();
972

973 974
        ParticleGroupIterator iterLeaves = _particleBlocks.begin();
        const ParticleGroupIterator iterEndLeaves = _particleBlocks.end();
975 976

        while(iterCells != iterEndCells && iterLeaves != iterEndLeaves){
977 978 979 980 981 982 983

            (*iterCells)->forEachCell(
                [&](SymbolCellClass* symb,
                    PoleCellClass* mult,
                    LocalCellClass* loc)
                {
                const int leafIdx = (*iterLeaves)->getLeafIndex(symb->getMortonIndex());
984 985
                FAssertLF(leafIdx != -1);
                ParticlesAttachedClass aLeaf = (*iterLeaves)->template getLeaf <ParticlesAttachedClass>(leafIdx);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
986
                FAssertLF(aLeaf.isAttachedToSomething());
987
                function(symb, mult, loc, &aLeaf);
988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003
            });

            ++iterCells;
            ++iterLeaves;
        }

        FAssertLF(iterCells == iterEndCells && iterLeaves == iterEndLeaves);
    }



    /** @brief, for statistic purpose, display each block with number of
   * cell, size of header, starting index, and ending index
   */
    void printInfoBlocks(){
        std::cout << "Group Tree information:\n";
1004 1005 1006
        std::cout << "\t Group Size = " << _nbElementsPerBlock << "\n";
        std::cout << "\t Tree height = " << _treeHeight << "\n";
        for(int idxLevel = 1 ; idxLevel < _treeHeight ; ++idxLevel){
1007
            std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
1008 1009
            std::cout << "Level " << idxLevel << ", there are " << levelBlocks.size() << " groups.\n";