FGroupTree.hpp 81.5 KB
Newer Older
1 2 3
// ==== CMAKE =====
//
// ================
4

5
//
6 7
#ifndef FGROUPTREE_HPP
#define FGROUPTREE_HPP
BRAMAS Berenger's avatar
BRAMAS Berenger committed
8
#include <vector>
9
#include <functional>
10

COULAUD Olivier's avatar
COULAUD Olivier committed
11 12 13 14 15
#include "Utils/FAssert.hpp"
#include "Utils/FPoint.hpp"
#include "Utils/FQuickSort.hpp"
#include "Containers/FTreeCoordinate.hpp"
#include "Containers/FCoordinateComputer.hpp"
16 17
#include "FGroupOfCells.hpp"
#include "FGroupOfParticles.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
18
#include "FGroupAttachedLeaf.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
19
#include "Kernels/P2P/FP2PParticleContainer.hpp"
20 21 22
#ifdef SCALFMM_USE_MPI
#include "FDistributedGroupTreeBuilder.hpp"
#endif
23

24
template <class FReal, class SymbolCellClass, class PoleCellClass, class LocalCellClass,
25
          class GroupAttachedLeafClass, unsigned NbSymbAttributes, unsigned NbAttributesPerParticle, class AttributeClass = FReal>
26
class FGroupTree {
27
public:
28 29 30 31 32 33 34
  typedef GroupAttachedLeafClass     BasicAttachedClass;   // Leaf
  typedef FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle,AttributeClass>     ParticleGroupClass;
  typedef FGroupOfCells<SymbolCellClass, PoleCellClass, LocalCellClass> CellGroupClass;
  typedef SymbolCellClass  GroupSymbolCellClass_T ;
  typedef LocalCellClass   GroupCellDownClass_T ;
  typedef PoleCellClass    GroupCellUpClass_T ;
  typedef GroupAttachedLeafClass     LeafClass_T;   // Leaf
BRAMAS Berenger's avatar
BRAMAS Berenger committed
35

36
protected:
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  //< height of the tree (1 => only the root)
  const int _treeHeight;
  //< max number of cells in a block
  const int _nbElementsPerBlock;
  //< all the blocks of the tree
  std::vector<CellGroupClass*>* _cellBlocksPerLevel;
  //< all the blocks of leaves
  std::vector<ParticleGroupClass*> _particleBlocks;

  //< the space system center
  const FPoint<FReal> boxCenter;
  //< the space system corner (used to compute morton index)
  const FPoint<FReal> boxCorner;
  //< the space system width
  const FReal boxWidth;
  //< the width of a box at width level
  const FReal boxWidthAtLeafLevel;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
54

55
public:
56 57 58 59
  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;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
60

61
  /** This constructor create a blocked octree from a usual octree
62 63 64 65
     * 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.
     */
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
  template<class OctreeClass>
  FGroupTree()
  {}
  template<class OctreeClass>
  FGroupTree(const int in_treeHeight, const int in_nbElementsPerBlock, OctreeClass*const inOctreeSrc)
    : _treeHeight(in_treeHeight), _nbElementsPerBlock(in_nbElementsPerBlock), _cellBlocksPerLevel(nullptr),
      boxCenter(inOctreeSrc->getBoxCenter()), boxCorner(inOctreeSrc->getBoxCenter(),-(inOctreeSrc->getBoxWidth()/2)),
      boxWidth(inOctreeSrc->getBoxWidth()), boxWidthAtLeafLevel(inOctreeSrc->getBoxWidth()/FReal(1<<(in_treeHeight-1))){

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

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

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

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
          // Create a block with the apropriate parameters
          CellGroupClass*const newBlock = new CellGroupClass(blockIteratorInOctree.getCurrentGlobalIndex(),
                                                             octreeIterator.getCurrentGlobalIndex()+1,
                                                             sizeOfBlock);
          FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<FReal, NbSymbAttributes, NbAttributesPerParticle, AttributeClass>(blockIteratorInOctree.getCurrentGlobalIndex(),
                                                                                                                                                                                                               octreeIterator.getCurrentGlobalIndex()+1,
                                                                                                                                                                                                               sizeOfBlock, nbParticlesInGroup);

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

              SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
              symbolic.setMortonIndex(newNodeIndex);
              symbolic.setCoordinate(newNodeCoordinate);
              symbolic.setLevel(idxLevel);

              // Add leaf
              nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(newNodeIndex, cellIdInBlock,
                                                                      blockIteratorInOctree.getCurrentLeaf()->getSrc()->getNbParticles(),
                                                                      nbParticlesOffsetBeforeLeaf);

              BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(cellIdInBlock);
              attachedLeaf.copyFromContainer(blockIteratorInOctree.getCurrentLeaf()->getSrc(), 0);

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

129 130 131
          // Keep the block
          _cellBlocksPerLevel[idxLevel].push_back(newBlock);
          _particleBlocks.push_back(newParticleBlock);
132

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

136 137
      avoidGotoLeft.moveUp();
      octreeIterator = avoidGotoLeft;
138 139
    }

140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
    // For each level from heigth - 2 to 1
    for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
        typename OctreeClass::Iterator avoidGotoLeft = octreeIterator;
        // For each cell at this level
        do {
            typename OctreeClass::Iterator blockIteratorInOctree = octreeIterator;
            // Move the iterator per _nbElementsPerBlock (or until it cannot move right)
            int sizeOfBlock = 1;
            while(sizeOfBlock < _nbElementsPerBlock && octreeIterator.moveRight()){
                sizeOfBlock += 1;
              }

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

                SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                symbolic.setMortonIndex(newNodeIndex);
                symbolic.setCoordinate(newNodeCoordinate);
                symbolic.setLevel(idxLevel);

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

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

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

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

  /**
184 185 186 187 188 189 190
     * 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
     */
191 192 193 194 195 196 197 198 199 200 201 202
  template<class ParticleContainer>
  FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
             const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
             const bool particlesAreSorted = false, MortonIndex inLeftLimite = -1):
    _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
    boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
    boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
  {
    _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];

    MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
    // First we work at leaf level
203
    {
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
      // 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,
                                                                                                        _treeHeight,
                                                                                                        FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
            const MortonIndex particleIndex = host.getMortonIndex();
            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;
            });
        }
233

234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
      FAssertLF(nbParticles == 0 || inLeftLimite < particlesToSort[0].mindex);
      // Convert to block
      const int idxLevel = (_treeHeight - 1);
      FSize* nbParticlesPerLeaf = new FSize[_nbElementsPerBlock];
      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
          while(sizeOfBlock < _nbElementsPerBlock && lastParticle < nbParticles){
              if(sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != particlesToSort[lastParticle].mindex){
                  currentBlockIndexes[sizeOfBlock] = particlesToSort[lastParticle].mindex;
                  nbParticlesPerLeaf[sizeOfBlock]  = 1;
                  sizeOfBlock += 1;
                }
              else{
                  nbParticlesPerLeaf[sizeOfBlock-1] += 1;
252
                }
253
              lastParticle += 1;
254
            }
255 256 257
          while(lastParticle < nbParticles && currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex){
              nbParticlesPerLeaf[sizeOfBlock-1] += 1;
              lastParticle += 1;
258 259
            }

260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
          // 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;
          FSize offsetParticles = firstParticle;
          for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
              newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

              SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
              symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
              FTreeCoordinate coord;
              coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
              symbolic.setCoordinate(coord);
              symbolic.setLevel(idxLevel);

              // 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);
289
                }
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
              offsetParticles += nbParticlesPerLeaf[cellIdInBlock];
            }

          // Keep the block
          _cellBlocksPerLevel[idxLevel].push_back(newBlock);
          _particleBlocks.push_back(newParticleBlock);

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

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

        CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
        const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();

        // 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
            while(sizeOfBlock < _nbElementsPerBlock && iterChildCells != iterChildEndCells ){
                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();
                      }
                  }
              }
347

348 349
            // If group is full
            if(sizeOfBlock == _nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
350 351
                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
352 353
                    currentBlockIndexes[sizeOfBlock-1]+1,
                    sizeOfBlock);
354 355 356 357
                // Init cells
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

358 359
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
360
                    FTreeCoordinate coord;
361
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
362 363
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
364
                  }
365 366

                // Keep the block
367
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
368 369

                sizeOfBlock = 0;
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 410 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 441 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 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498
              }
          }
      }
    delete[] currentBlockIndexes;
  }

  /**
     * 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>
  FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
             const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
             const bool particlesAreSorted, const bool oneParent,
             const FReal inCoverRatio = 0.0, MortonIndex inLeftLimite = -1):
    _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
    boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
    boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
  {

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

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

    MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
    // 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,
                                                                                                        _treeHeight,
                                                                                                        FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
            const MortonIndex particleIndex = host.getMortonIndex();
            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;
            });
        }

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

      // Convert to block
      const int idxLevel = (_treeHeight - 1);
      int* nbParticlesPerLeaf = new int[_nbElementsPerBlock];
      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
          while(sizeOfBlock < _nbElementsPerBlock && lastParticle < nbParticles
                && (userCoverRatio == false
                    || sizeOfBlock == 0
                    || currentBlockIndexes[sizeOfBlock-1] == particlesToSort[lastParticle].mindex
                    || (FReal(sizeOfBlock+1)/FReal(particlesToSort[lastParticle].mindex-particlesToSort[firstParticle].mindex)) >= inCoverRatio)){
              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);

              SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
              symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
              FTreeCoordinate coord;
              coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
              symbolic.setCoordinate(coord);
              symbolic.setLevel(idxLevel);

              // 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];
499
            }
500 501 502 503 504 505 506

          // Keep the block
          _cellBlocksPerLevel[idxLevel].push_back(newBlock);
          _particleBlocks.push_back(newParticleBlock);

          sizeOfBlock = 0;
          firstParticle = lastParticle;
507
        }
508 509 510
      delete[] nbParticlesPerLeaf;
      delete[] particlesToSort;
    }
511 512


513 514 515
    // For each level from heigth - 2 to 1
    for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
        inLeftLimite = (inLeftLimite == -1 ? inLeftLimite : (inLeftLimite>>3));
516

517 518
        CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
        const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
519

520 521 522 523 524 525 526 527 528 529 530 531 532
        // 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;
533

534
        if(oneParent == false){
535 536 537
            // 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
538
                while(sizeOfBlock < _nbElementsPerBlock && iterChildCells != iterChildEndCells ){
539
                    if((sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != (currentCellIndex>>3))
540
                       && (*iterChildCells)->exists(currentCellIndex)){
541 542 543
                        currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                        sizeOfBlock += 1;
                        currentCellIndex = (((currentCellIndex>>3)+1)<<3);
544
                      }
545 546
                    else{
                        currentCellIndex += 1;
547
                      }
548 549 550 551 552 553
                    // 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();
554 555 556
                          }
                      }
                  }
557 558

                // If group is full
559
                if(sizeOfBlock == _nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
560 561
                    // Create a group
                    CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
562 563
                        currentBlockIndexes[sizeOfBlock-1]+1,
                        sizeOfBlock);
564 565 566 567
                    // Init cells
                    for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                        newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

568 569
                        SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                        symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
570
                        FTreeCoordinate coord;
571
                        coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
572 573
                        symbolic.setCoordinate(coord);
                        symbolic.setLevel(idxLevel);
574
                      }
575 576

                    // Keep the block
577
                    _cellBlocksPerLevel[idxLevel].push_back(newBlock);
578 579

                    sizeOfBlock = 0;
580 581 582 583
                  }
              }
          }
        else{
584
            // We need to proceed each group in sub level
585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
            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)
                    while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
                        ++iterChildCells;
                        // Update morton index
                        if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                            currentCellIndex = (*iterChildCells)->getStartingIndex();
                          }
                      }
                  }
613

614 615 616 617
                // If group is full
                if(sizeOfBlock){
                    // Create a group
                    CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
618
                        currentBlockIndexes[sizeOfBlock-1]+1,
619 620 621 622
                        sizeOfBlock);
                    // Init cells
                    for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                        newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);
623

624 625 626 627 628 629 630
                        SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                        symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
                        FTreeCoordinate coord;
                        coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
                        symbolic.setCoordinate(coord);
                        symbolic.setLevel(idxLevel);
                      }
631

632 633
                    // Keep the block
                    _cellBlocksPerLevel[idxLevel].push_back(newBlock);
634

635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
                    sizeOfBlock = 0;
                  }
              }
          }
      }
    delete[] currentBlockIndexes;
  }
    /**
     * Sequential Constructor of GroupTree
     * used to construct a duplicated Ggroup tree on all processes
     * @param[in] in_treeHeight size of the tree
     * @param[in] in_boxWidth   bow witdh
     * @param[in] in_boxCenter  box center
     * @param[in] in__nbElementsPerBlock block size
     * @param[in] inParticlesContainer  an array of particles
     * @param[out] blockSizeAtEachLevel  box width at leaf level
     * @param[in] particlesAreSorted  True if the particle are sorted
     */
  template<class ParticleContainer>
  FGroupTree(const int in_treeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
             const int in_nbElementsPerBlock, ParticleContainer* inParticlesContainer,
             std::vector<std::vector<int>> & blockSizeAtEachLevel,
             const bool particlesAreSorted = false):
    _treeHeight(in_treeHeight),_nbElementsPerBlock(in_nbElementsPerBlock),_cellBlocksPerLevel(nullptr),
    boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
    boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(in_treeHeight-1)))
  {
    _cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];

    MortonIndex* currentBlockIndexes = new MortonIndex[_nbElementsPerBlock];
    // 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,
                                                                                                        _treeHeight,
                                                                                                        FPoint<FReal>(xpos[idxPart], ypos[idxPart], zpos[idxPart]) );
684 685
        //    const MortonIndex particleIndex         = host.getMortonIndex();
            particlesToSort[idxPart].mindex         = host.getMortonIndex();
686 687 688 689
            particlesToSort[idxPart].originalIndex = idxPart;
          }
      }

690

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

698 699 700 701 702 703 704 705 706 707 708 709 710 711 712
      // Convert to block
      const int idxLevel = (_treeHeight - 1);
      int idxBlock = 0;
      FSize* nbParticlesPerLeaf = new FSize[_nbElementsPerBlock];
      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
          while(sizeOfBlock < blockSizeAtEachLevel[_treeHeight-1][idxBlock] && lastParticle < nbParticles){
              if(sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != particlesToSort[lastParticle].mindex){
                  currentBlockIndexes[sizeOfBlock] = particlesToSort[lastParticle].mindex;
                  nbParticlesPerLeaf[sizeOfBlock]  = 1;
                  sizeOfBlock += 1;
713
                }
714 715
              else{
                  nbParticlesPerLeaf[sizeOfBlock-1] += 1;
716
                }
717
              lastParticle += 1;
718
            }
719 720 721 722 723 724 725 726 727 728 729 730 731 732
          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);

          /////////////////////////  TO REMOVE ?? //////////////
733 734 735 736
//          #include <iostream>
//	  using namespace std;
//	  if(currentBlockIndexes[sizeOfBlock-1]+1 == 511)
//	    cout << "Suricate" << endl;
737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
	  /////////////////////////////////////////////////////

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

	      SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
	      symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
	      FTreeCoordinate coord;
	      coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
	      symbolic.setCoordinate(coord);
	      symbolic.setLevel(idxLevel);

	      // 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
765
	  newBlock->declare_mine();
766 767 768 769 770 771 772 773 774
	  _cellBlocksPerLevel[idxLevel].push_back(newBlock);
	  _particleBlocks.push_back(newParticleBlock);

	  sizeOfBlock = 0;
	  firstParticle = lastParticle;
	  ++idxBlock;
	}
      delete[] nbParticlesPerLeaf;
      delete[] particlesToSort;
775
    }
776
//    MPI_Barrier(MPI_COMM_WORLD);
777

778 779 780 781
    // For each level from heigth - 2 to 1
    for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
        CellGroupConstIterator iterChildCells = _cellBlocksPerLevel[idxLevel+1].begin();
        const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
782

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

788 789 790
        MortonIndex currentCellIndex = (*iterChildCells)->getStartingIndex();
        int sizeOfBlock = 0;
        int idxBlock    = 0;
791

792 793
        // We need to proceed each group in sub level
        while(iterChildCells != iterChildEndCells){
794

795 796
            // Count until end of sub group is reached or we have enough cells
            while(sizeOfBlock < blockSizeAtEachLevel[idxLevel][idxBlock] && iterChildCells != iterChildEndCells ){
797

798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815
                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();
                      }
                  }
              }
816

817 818
            // If group is full
            if(sizeOfBlock == blockSizeAtEachLevel[idxLevel][idxBlock] || (sizeOfBlock && iterChildCells == iterChildEndCells)){ //NOTE la seconde partie va sûrement sauter, car la taille est pré-calculée
819 820
                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
821 822
                    currentBlockIndexes[sizeOfBlock-1]+1,
                    sizeOfBlock);
823 824 825 826
                // Init cells
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);

827 828
                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
829
                    FTreeCoordinate coord;
830 831 832
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
833
                  }
834 835

                // Keep the block
836
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
837 838

                sizeOfBlock = 0;
839 840 841 842 843 844
                ++idxBlock;
              }
          }
      }
    delete[] currentBlockIndexes;
  }
845 846


847
  /**
848 849 850 851 852 853 854 855 856
     * 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
     */
857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877
  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)
  {
    this->_cellBlocksPerLevel = new std::vector<CellGroupClass*>[_treeHeight];

  }


  /**
878 879 880 881 882 883
     * 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
     */
884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
  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
  /////////////////////////////////////////////////////////

  /**
916 917 918 919 920 921 922
     * 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
     */
923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941
  template<class Group_Linear_tree,
           class Particle_Container
           >
  void create_tree(Group_Linear_tree& in_lin_tree,
                   const Particle_Container& particles){
    MortonIndex in_left_limit = in_lin_tree.get_left_limit();
    // Creation of the leaf level and groups of particle
    std::vector<MortonIndex> current_block_indexes = create_leaf_level(in_lin_tree,particles);
    // Creation of every level of the tree
    create_block_nodes_level(
          current_block_indexes,
          in_left_limit);
  }

  /** This function dealloc the tree by deleting each block */
  ~FGroupTree(){
    for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        for (CellGroupClass* block: levelBlocks){
942
            delete block;
943 944 945 946 947 948 949 950 951 952 953 954 955
          }
      }
    delete[] _cellBlocksPerLevel;
    for (ParticleGroupClass* block: _particleBlocks){
        delete block;
      }
  }

  ////////////////////////////////////////////////////////
  // Lambda function to apply to all member
  /////////////////////////////////////////////////////////

  /**
956 957 958
   * @brief forEachLeaf iterate on the leaf and apply the function
   * @param function
   */
959 960 961 962 963 964 965 966
  template<class ParticlesAttachedClass>
  void forEachLeaf(std::function<void(ParticlesAttachedClass*)> function){
    for (ParticleGroupClass* block: _particleBlocks){
        block->forEachLeaf(function);
      }
  }
  /**
   * @brief forEachMyLeaf iterate on the leaf and apply the function
967 968
   * @param function
   */
969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
  template<class ParticlesAttachedClass>
  void forEachMyLeaf(std::function<void(ParticlesAttachedClass*)> function){
    for (ParticleGroupClass* block: _particleBlocks){
        if(block->isMine())
          block->forEachLeaf(function);
      }
  }

  /**
   * @brief forEachCell iterate on the cell and apply the function
   * @param function
   */
  void forEachCell(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*)> function){
    for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        for (CellGroupClass* block: levelBlocks){
            block->forEachCell(function);
          }
      }
  }

  void forEachMyCell(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*)> function){
    for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        for (CellGroupClass* block: levelBlocks){
            if(block->isMine())
              block->forEachCell(function);
          }
      }
  }

  /**
1001 1002 1003
   * @brief forEachLeaf iterate on the cell and apply the function
   * @param function
   */
1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023
  void forEachCellWithLevel(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,const int)> function){
    for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        for (CellGroupClass* block: levelBlocks){
            block->forEachCell(function, idxLevel);
          }
      }
  }

  void forEachMyCellWithLevel(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,const int)> function){
    for(int idxLevel = 0 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        for (CellGroupClass* block: levelBlocks){
            if(block->isMine())
              block->forEachCell(function, idxLevel);
          }
      }
  }

  /**
1024 1025 1026 1027 1028 1029 1030
   * @brief forEachLeaf iterate on the cells that i own and apply the function on the leaves
   *
   *  We iterate on the owned group and for for each group we iterate on its cells.
   * We obtain the leaf through the morton index of the cell and we apply the function on it.
   *
   * @param function function to apply on each leaf
   *
1031
   */
1032 1033
  template<class ParticlesAttachedClass>
  void forEachCellLeaf(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,ParticlesAttachedClass*)> function){
1034
    CellGroupIterator       iterCells    = _cellBlocksPerLevel[_treeHeight-1].begin();
1035
    const CellGroupIterator iterEndCells = _cellBlocksPerLevel[_treeHeight-1].end();
1036

1037
    ParticleGroupIterator       iterLeaves    = _particleBlocks.begin();
1038
    const ParticleGroupIterator iterEndLeaves = _particleBlocks.end();
1039 1040 1041 1042 1043
   // int count = 0 ;
     //  Iterate on Cell group and Leaf group
   while(iterCells != iterEndCells && iterLeaves != iterEndLeaves){
       //
       //  Iterate all Cells inside the group
1044 1045 1046 1047 1048
        (*iterCells)->forEachCell(
              [&](SymbolCellClass* symb,
              PoleCellClass* mult,
              LocalCellClass* loc)
        {
1049 1050 1051 1052 1053
         //  The group of cells is mine --> Leaf exist
//        std::cout << " forEachCellLeaf " << (*iterCells)->isMine() << " morton " << symb->getMortonIndex()
//                  << "   leafIdx " << (*iterLeaves)->getLeafIndex(symb->getMortonIndex())
//                  << "   count  " <<count <<std::endl;
        if ((*iterCells)->isMine() ) {
1054

1055
              //  get leafindex (position in the group) from the Morton index
1056
              const int leafIdx = (*iterLeaves)->getLeafIndex(symb->getMortonIndex());
1057
              //  Leaf exists and we apply function on it
1058 1059 1060 1061
              FAssertLF(leafIdx != -1);
              ParticlesAttachedClass aLeaf = (*iterLeaves)->template getLeaf <ParticlesAttachedClass>(leafIdx);
              FAssertLF(aLeaf.isAttachedToSomething());
              function(symb, mult, loc, &aLeaf);
1062 1063 1064
            }
   //       ++count;
        });
1065 1066 1067
        ++iterCells;
        ++iterLeaves;
      }
1068

1069 1070
    FAssertLF(iterCells == iterEndCells && iterLeaves == iterEndLeaves);
  }
1071 1072


1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103
//  template<class ParticlesAttachedClass>
//  void forEachCellMyLeaf(std::function<void(SymbolCellClass*,PoleCellClass*,LocalCellClass*,ParticlesAttachedClass*)> function){
//    CellGroupIterator iterCells = _cellBlocksPerLevel[_treeHeight-1].begin();

//    const CellGroupIterator iterEndCells = _cellBlocksPerLevel[_treeHeight-1].end();

//    ParticleGroupIterator iterLeaves = _particleBlocks.begin();
//    const ParticleGroupIterator iterEndLeaves = _particleBlocks.end();

//    while(iterCells != iterEndCells && iterLeaves != iterEndLeaves){
//        if((*iterCells)->isMine()){
//            (*iterCells)->forEachCell(
//                  [&](SymbolCellClass* symb,
//                  PoleCellClass*       mult,
//                  LocalCellClass*      loc)
//            {
//              const int leafIdx = (*iterLeaves)->getLeafIndex(symb->getMortonIndex());
//              FAssertLF(leafIdx != -1);
//              ParticlesAttachedClass aLeaf = (*iterLeaves)->template getLeaf <ParticlesAttachedClass>(leafIdx);
//              FAssertLF(aLeaf.isAttachedToSomething());
//              function(symb, mult, loc, &aLeaf);
//            });
//          }
//        ++iterCells;
//        ++iterLeaves;
//      }

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


1104
  /** @brief, for statistic purpose, display each block with number of
1105 1106
   * cell, size of header, starting index, and ending index
   */
1107 1108 1109 1110 1111 1112 1113
  void printInfoBlocks(){
    std::cout << "Group Tree information:\n";
    std::cout << "\t Group Size = " << _nbElementsPerBlock << "\n";
    std::cout << "\t Tree height = " << _treeHeight << "\n";
    for(int idxLevel = 1 ; idxLevel < _treeHeight ; ++idxLevel){
        std::vector<CellGroupClass*>& levelBlocks = _cellBlocksPerLevel[idxLevel];
        std::cout << "Level " << idxLevel << ", there are " << levelBlocks.size() << " groups.\n";
1114
        int idxGroup = 0;
1115
        for (const CellGroupClass* block: levelBlocks){
1116
            std::cout << "\t Group " << (idxGroup++);
1117 1118
          //  std::cout << "\t local " << std::boolalpha << block->isMine();
            std::cout << "\t Size = " << block->getNumberOfCellsInBlock();
1119 1120
            std::cout << "\t Starting Index = " << block->getStartingIndex();
            std::cout << "\t Ending Index = " << block->getEndingIndex();
1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238
         //   std::cout << "\t Global index  = " << block->getIdxGlobal();
            std::cout << "\t Ratio of usage = " <<
                         float(block->getNumberOfCellsInBlock())/float(block->getEndingIndex()-block->getStartingIndex()) << "\n";
          }
      }

    std::cout << "There are " << _particleBlocks.size() << " leaf-groups.\n";
    int idxGroup = 0;
    FSize totalNbParticles = 0;
    for (const ParticleGroupClass* block: _particleBlocks){
        std::cout << "\t Group " << (idxGroup++);

        std::cout << "\t Size = " << block->getNumberOfLeavesInBlock();
        std::cout << "\t Starting Index = " << block->getStartingIndex();
        std::cout << "\t Ending Index = " << block->getEndingIndex();
        std::cout << "\t Nb Particles = " << block->getNbParticlesInGroup();
        std::cout << "\t Global index  = " << block->getIdxGlobal();
        std::cout << "\t Ratio of usage = " << float(block->getNumberOfLeavesInBlock())/float(block->getEndingIndex()-block->getStartingIndex()) << "\n";
        totalNbParticles += block->getNbParticlesInGroup();
      }
    std::cout << "There are " << totalNbParticles << " particles.\n";
  }

  /////////////////////////////////////////////////////////
  // Algorithm function
  /////////////////////////////////////////////////////////

  int getHeight() const {
    return _treeHeight;
  }

  CellGroupIterator cellsBegin(const int inLevel){
    FAssertLF(inLevel < _treeHeight);
    return _cellBlocksPerLevel[inLevel].begin();
  }

  CellGroupConstIterator cellsBegin(const int inLevel) const {
    FAssertLF(inLevel < _treeHeight);
    return _cellBlocksPerLevel[inLevel].begin();
  }

  CellGroupIterator cellsEnd(const int inLevel){
    FAssertLF(inLevel < _treeHeight);
    return _cellBlocksPerLevel[inLevel].end();
  }

  CellGroupConstIterator cellsEnd(const int inLevel) const {
    FAssertLF(inLevel < _treeHeight);
    return _cellBlocksPerLevel[inLevel].end();
  }

  int getNbCellGroupAtLevel(const int inLevel) const {
    FAssertLF(inLevel < _treeHeight);
    return int(_cellBlocksPerLevel[inLevel].size());
  }

  CellGroupClass* getCellGroup(const int inLevel, const int inIdx){
    FAssertLF(inLevel < _treeHeight);
    FAssertLF(inIdx < int(_cellBlocksPerLevel[inLevel].size()));
    return _cellBlocksPerLevel[inLevel][inIdx];
  }

  const int getNbElementsPerBlock() const{
    return this->_nbElementsPerBlock;
  }

  const CellGroupClass* getCellGroup(const int inLevel, const int inIdx) const {
    FAssertLF(inLevel < _treeHeight);
    FAssertLF(inIdx < int(_cellBlocksPerLevel[inLevel].size()));
    return _cellBlocksPerLevel[inLevel][inIdx];
  }

  ParticleGroupIterator leavesBegin(){
    return _particleBlocks.begin();
  }

  ParticleGroupConstIterator leavesBegin() const {
    return _particleBlocks.begin();
  }

  ParticleGroupIterator leavesEnd(){
    return _particleBlocks.end();
  }

  ParticleGroupConstIterator leavesEnd() const {
    return _particleBlocks.end();
  }

  int getNbParticleGroup() const {
    return int(_particleBlocks.size());
  }

  ParticleGroupClass* getParticleGroup(const int inIdx){
    FAssertLF(inIdx < int(_particleBlocks.size()));
    return _particleBlocks[inIdx];
  }

  const ParticleGroupClass* getParticleGroup(const int inIdx) const {
    FAssertLF(inIdx < int(_particleBlocks.size()));
    return _particleBlocks[inIdx];
  }

  const FPoint<FReal> getBoxCenter() const{
    return this->boxCenter;
  }

  const FReal getBoxWidth() const{
    return this->boxWidth;
  }

  std::size_t getTotalNbLeaf() {
    std::size_t nbLeaf = 0;
    for(int i = 0 ; i < this->getNbParticleGroup();++i){
        nbLeaf += this->_particleBlocks[i]->getNumberOfLeavesInBlock();
      }
    return nbLeaf;
  }
  /**
1239 1240 1241 1242 1243 1244 1245
     * RESTRICTION : The array will be initialise BEFORE
     * RESTRICTION : The morton index of particle will be at _treeHeight
     * get_number_of_particle compute the total number of
     * particle on every leaf, he just fill the array nb_particles_per_leaf
     * @author benjamin.dufoyer@inria.fr
     * @param  container    container of particle
     */
1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265
  template<class particle_t>
  void get_number_of_particle(const std::vector<particle_t>& container,
                              std::vector<std::size_t>& nb_particles_per_leaf){
    FAssert(container.size() != 0);
    int current_idx = 0;
    std::size_t old_m_index   = container.front().morton_index;
    std::size_t current_m_idx = old_m_index;
    for(std::size_t i = 0 ; i < container.size(); ++i){
        current_m_idx = container[i].morton_index;
        if(current_m_idx == old_m_index){
            nb_particles_per_leaf[current_idx] += 1;
          } else {
            current_idx += 1;
            nb_particles_per_leaf[current_idx] += 1;
            old_m_index = current_m_idx;
          }
      }
  }

  /**
1266 1267 1268 1269 1270 1271 1272
     * create_leaf_level create the leaf level of the
     * Group tree from a blocked linear tree
     * @author benjamin.dufoyer@inria.fr
     * @param  in_lin_tree  Blocked linear tree
     * @param  particles    container of particle, will be a std::vector
     */

1273 1274 1275 1276 1277 1278
  template<class Blocked_Linear_tree,
           class Particle_Container>
  std::vector<MortonIndex> create_leaf_level(Blocked_Linear_tree& in_lin_tree,
                                             Particle_Container& particles)
  {
    // set parametter for the function
1279
    const int idxLevel = this->_treeHeight-1;
1280
    const int nb_block = in_lin_tree.get_nb_block();
1281
    const int block_size = this->_nbElementsPerBlock;
1282 1283 1284 1285 1286 1287 1288 1289
    std::size_t in_nb_leaf = in_lin_tree.get_nb_leaf();
    auto tree = in_lin_tree.get_tree();
    // alloc the vector for the current block index
    // get the number of particle per leaf
    std::vector<MortonIndex> current_block_indexes(this->_nbElementsPerBlock,0);
    std::vector<std::size_t> nb_particle_per_leaf(in_nb_leaf,0);
    this->get_number_of_particle(particles,nb_particle_per_leaf);
    // put the particle in the FP2PParticleContainer
1290
    FP2PParticleContainer<FReal> particle_container;
1291 1292 1293
    for(unsigned i = 0 ; i < particles.size() ; ++i){
        particle_container.push(particles[i].position(), particles[i].physicalValue());
      }
1294

1295 1296
    std::size_t leaf_number = 0;
    std::size_t leaf_number_min = 0;
1297 1298

    // Create every block
1299
    std::size_t idx_particules = 0;
1300 1301 1302 1303
    for(int n_block = 0 ; n_block < nb_block ; ++n_block){
        // Compute the morton index for the first and the
        // last cell of the block
        unsigned size_of_block = 0;
1304 1305 1306
        while(size_of_block < (unsigned)block_size
              && leaf_number < in_nb_leaf)
          {
1307 1308 1309
            current_block_indexes[size_of_block] = tree->data()[leaf_number].morton_index;
            leaf_number += 1;
            size_of_block += 1;
1310
          }
1311 1312

        CellGroupClass*const new_block = new CellGroupClass(current_block_indexes[0],
1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405
            current_block_indexes[size_of_block-1]+1, //+1 is need by the class
            size_of_block);
        size_t current_nb_particle = 0;
        for(size_t i = 0 ; i < size_of_block ; ++i){
            current_nb_particle += nb_particle_per_leaf[leaf_number_min+i];
          }
        FGroupOfParticles<
            FReal,
            NbSymbAttributes,
            NbAttributesPerParticle,
            AttributeClass>*const new_particle_block
            = new FGroupOfParticles<
            FReal,
            NbSymbAttributes,
            NbAttributesPerParticle,
            AttributeClass>
            (current_block_indexes[0],
            current_block_indexes[size_of_block-1]+1,
            size_of_block,
            current_nb_particle);

        // Initialise each cell of the block
        size_t nb_particles_offset_before_leaf = 0;
        for(unsigned cell_id_in_block = 0;  cell_id_in_block < size_of_block; ++cell_id_in_block)
          {
            // Adding cell into leaf block
            new_block->newCell(
                  current_block_indexes[cell_id_in_block],
                  cell_id_in_block);

            // Fill symbolic information of the block
            SymbolCellClass& symbolic =
                new_block->getSymbolic(cell_id_in_block);
            symbolic.setMortonIndex(current_block_indexes[cell_id_in_block]);
            FTreeCoordinate coord;
            coord.setPositionFromMorton(current_block_indexes[cell_id_in_block]);
            symbolic.setCoordinate(coord);
            symbolic.setLevel(idxLevel);

            // Adding cell into particle blockCells

            nb_particles_offset_before_leaf =
                new_particle_block->newLeaf(
                  current_block_indexes[cell_id_in_block],
                  cell_id_in_block,
                  FSize(nb_particle_per_leaf[leaf_number_min+cell_id_in_block]),
                nb_particles_offset_before_leaf
                );


            BasicAttachedClass attached_leaf =
                new_particle_block->template getLeaf<BasicAttachedClass>(cell_id_in_block);

            // Adding particle
            for(size_t idxPart = 0 ; idxPart <   nb_particle_per_leaf[leaf_number_min+cell_id_in_block] ; ++idxPart ){
                attached_leaf.setParticle(
                      idxPart,
                      idx_particules,
                      //nb_particles_offset_before_leaf+idxPart,
                      &particle_container);
                ++idx_particules;
              }
            // Setting the offset to don't use particle twice
            //offset_particles += nb_particle_per_leaf[idx_nb_particle_in_block];
            //idx_nb_particle_in_block += 1;
            // cell_id_in_block += 1;
          }
        leaf_number_min = leaf_number;
        new_block->declare_mine();
        //Stock the block cell and the block particles
        _cellBlocksPerLevel[idxLevel].push_back(new_block);
        _particleBlocks.push_back(new_particle_block);
        size_of_block = 0;
      }
    return {current_block_indexes.begin(),current_block_indexes.end()};

  }

  /**
    * create_level create every level
    * It's juste a factorisation from the Beregenger constructor
    * @author benjamin.dufoyer@inria.fr
    * @param currentBlockIndexes block repartition at leaf level
    * to construct
    * @param inLeftLimite left limit of block of the current proc
    * this parameter is not used with the blocked_linear_tree, he is here
    * to have compatibility with old constructor
    */
  void create_block_nodes_level(std::vector<MortonIndex>& currentBlockIndexes,
                                MortonIndex inLeftLimite = -1
      ){
    // Cronstruct every level
    for(int idxLevel = _treeHeight-2; idxLevel > 0 ; --idxLevel){
1406
        inLeftLimite = (inLeftLimite == -1 ? inLeftLimite : (inLeftLimite>>3));
1407

1408
        CellGroupConstIterator iterChildCells          = _cellBlocksPerLevel[idxLevel+1].begin();
1409
        const CellGroupConstIterator iterChildEndCells = _cellBlocksPerLevel[idxLevel+1].end();
1410 1411

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

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

        // We need to proceed each group in sub level
1426
        while(iterChildCells != iterChildEndCells){
1427
            // Count until end of sub group is reached or we have enough cells
1428 1429
            while(sizeOfBlock < _nbElementsPerBlock && iterChildCells != iterChildEndCells ){
                if((sizeOfBlock == 0 || currentBlockIndexes[sizeOfBlock-1] != (currentCellIndex>>3))
1430
                   && (*iterChildCells)->exists(currentCellIndex)){
1431 1432 1433
                    currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
                    sizeOfBlock += 1;
                    currentCellIndex = (((currentCellIndex>>3)+1)<<3);
1434
                  }
1435 1436
                else{
                    currentCellIndex += 1;
1437
                  }
1438
                // If we are at the end of the sub group, move to next
1439 1440
                while(iterChildCells != iterChildEndCells && (*iterChildCells)->getEndingIndex() <= currentCellIndex){
                    ++iterChildCells;
1441
                    // Update morton index
1442 1443
                    if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
                        currentCellIndex = (*iterChildCells)->getStartingIndex();
1444 1445 1446
                      }
                  }
              }
1447 1448

            // If group is full
1449 1450 1451
            if(sizeOfBlock == _nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
                // Create a group
                CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
1452 1453
                    currentBlockIndexes[sizeOfBlock-1]+1,
                    sizeOfBlock);
1454 1455 1456
                // Init cells
                for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
                    newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);
1457 1458

                    SymbolCellClass& symbolic = newBlock->getSymbolic(cellIdInBlock);
1459
                    symbolic.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
1460
                    FTreeCoordinate coord;
1461
                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock]);
1462 1463
                    symbolic.setCoordinate(coord);
                    symbolic.setLevel(idxLevel);
1464 1465
                  }
                newBlock->declare_mine();
1466
                // Keep the block
1467
                _cellBlocksPerLevel[idxLevel].push_back(newBlock);
1468
                sizeOfBlock = 0;
1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566
              }
          }
      }
  }



  /**
    * This function add all LET block put in parameter
    * She put block in order according to idx_global
    * She detect if we are at leaf level and create particle group
    *
    * @author benjamin.dufoyer@inria.fr
    * @param  block_to_insert pair of symbolic information of cellGroup and Part
    *                         group
    * @param  level           The level where are adding LET group
    */
  template<class particle_symbolic_block_t,
           class cell_symbolic_block_t>
  void add_LET_block(
      std::pair<std::vector<cell_symbolic_block_t>,
      std::vector<particle_symbolic_block_t>>& block_to_insert,
      int                                                level
      ){
    // Check if we are at leaf level
    bool leaf_level = ( level == ( _treeHeight - 1 ) );
    // Bind the vector of the pair
    std::vector<cell_symbolic_block_t> cell_to_insert = block_to_insert.first;
    std::vector<particle_symbolic_block_t> particle_to_insert = block_to_insert.second;
    // If we are at leaf level
    if(leaf_level){
        // Check if we have the same number of symoblic information of cellBlock
        // and of particleBlock
        FAssert(cell_to_insert.size() == particle_to_insert.size());
      } else {
        // Else check if the particle block is empty
        FAssert(particle_to_insert.size() == 0);
      }
    // if we have no block to insert, we don't need to continue this function
    if(cell_to_insert.size() == 0)
      return;
    // Get my local minimum index global
    int min_idx_global = this->getCellGroup(level,0)->getIdxGlobal();
    // Allocate vector of new block
    std::vector<CellGroupClass*> vect_block(cell_to_insert.size());

    // Fill the vector of new block
    unsigned block_at_begin = 0;
    // iterate on every cell
    for(unsigned i = 0; i < cell_to_insert.size(); ++i){
        // create new cell
        vect_block[i] = new CellGroupClass(
              cell_to_insert[i].start_index ,
              cell_to_insert[i].end_index,
              (int)cell_to_insert[i].nb_leaf_in_block );
        // set the global index of the cell
        vect_block[i]->setIdxGlobal(cell_to_insert[i].idx_global_block);
        // if the global index is less than the local idex, we need to
        // insert
        // the block at the beginning of the tree
        if(cell_to_insert[i].idx_global_block < min_idx_global){
            ++block_at_begin;
          }
        // init each cell of the new block
        for(unsigned j = 0; j < cell_to_insert[i].m_idx_in_block.size(); ++j){
            vect_block[i]->newCell(cell_to_insert[i].m_idx_in_block[j],j);
          }
      }
    // Add block at beginning of the level
    _cellBlocksPerLevel[level].insert(
          _cellBlocksPerLevel[level].begin(),
          vect_block.begin(),
          vect_block.begin()+block_at_begin);

    // Add block a the end of the level
    _cellBlocksPerLevel[level].insert(
          _cellBlocksPerLevel[level].end(),
          vect_block.begin()+block_at_begin,
          vect_block.end());
    // if we are at the leaf level
    if(leaf_level ){
        // init of the vector of particle
        std::vector<ParticleGroupClass*> vect_particle(particle_to_insert.size());

        block_at_begin = 0;
        // iterate on every symbolic particle group
        for(unsigned i = 0 ; i < particle_to_insert.size(); ++i ){
            // create a new particle group
            vect_particle[i] = new ParticleGroupClass(
                  cell_to_insert[i].start_index ,
                  cell_to_insert[i].end_index,
                  (int)cell_to_insert[i].nb_leaf_in_block,
                  particle_to_insert[i].nb_particles);
            // set the global index of the new particle group
            vect_particle[i]-> setIdxGlobal(particle_to_insert[i].idx_global_block);
            // if the current idx global block have a idx global smaller than
            // the global index in local
            if(cell_to_insert[i].idx_global_block < min_idx_global){
1567
                ++block_at_begin;
1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578
              }
            size_t offset = 0;
            // init all leaf of the current particle group
            for(int j = 0; j < cell_to_insert[i].nb_leaf_in_block; ++j){
                offset = vect_particle[i]->newLeaf(
                      cell_to_insert[i].m_idx_in_block[j],
                      j,
                      particle_to_insert[i].nb_particle_per_leaf[j],
                      offset);
              }
          }
1579
        // Add block at beginning of the level
1580 1581 1582 1583
        _particleBlocks.insert(
              _particleBlocks.begin(),
              vect_particle.begin(),
              vect_particle.begin()+block_at_begin);
DUFOYER Benjamin's avatar
DUFOYER Benjamin committed
1584
        // Add block a the end of the level
1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 1866
        _particleBlocks.insert(
              _particleBlocks.end(),
              vect_particle.begin()+block_at_begin,
              vect_particle.end());
      }
  }


#ifdef SCALFMM_USE_MPI
  /**
     * This function compute and add the local essential tree (LET) at
     * the level.
     * We compute interaction for the P2P(if needed) and M2L. We communicate
     * other proc to get the GroupOfCell needed for building the LET
     * @author benjamin.dufoyer@inria.fr
     * @param  group_linear_tree        The group linear tree
     * @param  level                    The level to build the LET
     * @param  dim                      The dimension of Coordinate
     */
  template<class GroupLinearTree>
  void create_LET_at_level(
      GroupLinearTree&    group_linear_tree,
      int&                level,
      MortonIndex&        gmin,
      MortonIndex&        gmax,
      MortonIndex&        lmin,
      MortonIndex&        lmax,
      int                 dim = 3
      ){
    // stock in the variable if we are at the leaf level
    bool leaf_level = (this->getHeight()-1 == level);
    // update the morton index
    if(!leaf_level){
        gmin = gmin >> 3;
        gmax = gmax >> 3;
      }
    const MortonIndex global_min_m_idx = gmin;
    const MortonIndex global_max_m_idx = gmax;
    // Compute min and max local morton index at the level needed
    if(this->getNbCellGroupAtLevel(level) > 0){
        lmin = this->getCellGroup(level,0)->getStartingIndex();
        lmax = this->getCellGroup(level,this->getNbCellGroupAtLevel(level)-1)->getEndingIndex()-1;
      } else {
        lmin = -1;
        lmax = -1;
      }
    const MortonIndex local_min_m_idx = lmin;
    const MortonIndex local_max_m_idx = lmax;

    // declare variable, needed because we fill it in a if case
    std::vector<MortonIndex> leaf_P2P;
    if(leaf_level){
        // IDEA : can be a task
        // This function compute the leaf needed by the P2P operation
        // This function return a vector with all leaf needed
        // The P2P interaction is only needed at leaf level
        leaf_P2P = dstr_grp_tree_builder::get_leaf_P2P_interaction(
              *this,
              global_min_m_idx,
              global_max_m_idx,
              local_min_m_idx,
              local_max_m_idx);
      }

    // IDEA can be a task
    // This function compute the leaf needed by the M2L operation
    // This function return a vector with all leaf needed
    // get leaf M2L
    std::vector<MortonIndex> leaf_M2L =
        dstr_grp_tre