Create real class for grouped tree

Src/GroupTree/FGroupAttachedLeaf.hpp

+#ifndef FGROUPATTACHEDLEAF_HPP
+#define FGROUPATTACHEDLEAF_HPP
+
+#include "../Utils/FGlobal.hpp"
+
+/**
+ * This class is "attached" to a buffer.
+ * It is a wrapper/interface to a memory area that stores all the particles info.
+ * The idea is to hidde the group allocation from the group tree but
+ * to keep the same interface than the FBasicParticlesContainer.
+ */
+template <unsigned NbAttributesPerParticle, class AttributeClass = FReal>
+class FGroupAttachedLeaf {
+protected:
+    //< Nb of particles in the current leaf
+    const int nbParticles;
+    //< Pointers to the positions of the particles
+    FReal* positionsPointers[3];
+    //< Pointers to the attributes of the particles
+    AttributeClass* attributes[NbAttributesPerParticle];
+
+    // Forbid copy even if there is no real reason to do that
+    FGroupAttachedLeaf(const FGroupAttachedLeaf&) = delete;
+    FGroupAttachedLeaf& operator=(const FGroupAttachedLeaf&) = delete;
+
+public:
+    /**
+     * @brief FGroupAttachedLeaf
+     * @param inNbParticles the number of particles in the leaf
+     * @param inPositionBuffer the memory address of the X array of particls
+     * @param inLeadingPosition each position is access by inPositionBuffer + in bytes inLeadingPosition*idx
+     * @param inAttributesBuffer the memory address of the first attribute
+     * @param inLeadingAttributes each attribute is access by inAttributesBuffer + in bytes inLeadingAttributes*idx
+     */
+    FGroupAttachedLeaf(const int inNbParticles, FReal* inPositionBuffer, const size_t inLeadingPosition,
+                       AttributeClass* inAttributesBuffer, const size_t inLeadingAttributes)
+        : nbParticles(inNbParticles){
+        // Redirect pointers to position
+        positionsPointers[0] = inPositionBuffer;
+        positionsPointers[1] = reinterpret_cast<FReal*>(reinterpret_cast<unsigned char*>(inPositionBuffer) + inLeadingPosition);
+        positionsPointers[2] = reinterpret_cast<FReal*>(reinterpret_cast<unsigned char*>(inPositionBuffer) + inLeadingPosition*2);
+
+        // Redirect pointers to data
+        for(int idxAttribute = 0 ; idxAttribute < NbAttributesPerParticle ; ++idxAttribute){
+            particleAttributes[idxAttribute] = reinterpret_cast<AttributeClass*>(reinterpret_cast<unsigned char*>(inAttributesBuffer) + idxAttribute*inLeadingAttributes);
+        }
+    }
+
+    /**
+     * @brief getNbParticles
+     * @return the number of particles in the leaf
+     */
+    int getNbParticles() const{
+        return nbParticles;
+    }
+
+    /**
+     * @brief getPositions
+     * @return a FReal*[3] to get access to the positions
+     */
+    const FReal*const* getPositions() const {
+        return positionsPointers;
+    }
+
+    /**
+     * @brief getWPositions
+     * @return get the position in write mode
+     */
+    FReal* const* getWPositions() {
+        return positionsPointers;
+    }
+
+    /**
+     * @brief getAttribute
+     * @param index
+     * @return the attribute at index index
+     */
+    AttributeClass* getAttribute(const int index) {
+        return attributes[index];
+    }
+
+    /**
+     * @brief getAttribute
+     * @param index
+     * @return
+     */
+    const AttributeClass* getAttribute(const int index) const {
+        return attributes[index];
+    }
+
+    /**
+     * Get the attribute with a forcing compile optimization
+     */
+    template <int index>
+    AttributeClass* getAttribute() {
+        static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
+        return attributes[index];
+    }
+
+    /**
+     * Get the attribute with a forcing compile optimization
+     */
+    template <int index>
+    const AttributeClass* getAttribute() const {
+        static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
+        return attributes[index];
+    }
+};
+
+#endif // FGROUPATTACHEDLEAF_HPP

Src/GroupTree/FGroupOfCells.hpp

+#ifndef FGROUPOFCELLS_HPP
+#define FGROUPOFCELLS_HPP
+
+#include "../Utils/FAssert.hpp"
+#include "../Containers/FTreeCoordinate.hpp"
+
+#include <list>
+#include <functional>
+
+
+/**
+* @brief The FGroupOfCells class manages the cells in block allocation.
+*/
+template <class CellClass>
+class FGroupOfCells {
+    /** One header is allocated at the beginning of each block */
+    struct BlockHeader{
+        MortonIndex startingIndex;
+        MortonIndex endingIndex;
+        int numberOfCellsInBlock;
+        int blockIndexesTableSize;
+    };
+
+protected:
+    //< Pointer to a block memory
+    unsigned char* memoryBuffer;
+
+    //< Pointer to the header inside the block memory
+    BlockHeader*    blockHeader;
+    //< Pointer to the indexes table inside the block memory
+    int*            blockIndexesTable;
+    //< Pointer to the cells inside the block memory
+    CellClass*      blockCells;
+    //< This value is for not used cells
+    static const MortonIndex CellIsEmptyFlag = -1;
+
+public:
+    /**
+ * @brief FGroupOfCells
+ * @param inStartingIndex first cell morton index
+ * @param inEndingIndex last cell morton index + 1
+ * @param inNumberOfCells total number of cells in the interval (should be <= inEndingIndex-inEndingIndex)
+ */
+    FGroupOfCells(const MortonIndex inStartingIndex, const MortonIndex inEndingIndex, const int inNumberOfCells)
+        : memoryBuffer(0), blockHeader(0), blockIndexesTable(0), blockCells(0) {
+        // Find the number of cell to allocate in the blocks
+        const int blockIndexesTableSize = int(inEndingIndex-inStartingIndex);
+        FAssertLF(inNumberOfCells <= blockIndexesTableSize);
+        // Total number of bytes in the block
+        const size_t memoryToAlloc = sizeof(BlockHeader) + (blockIndexesTableSize*sizeof(int)) + (inNumberOfCells*sizeof(CellClass));
+
+        // Allocate
+        memoryBuffer = new unsigned char[memoryToAlloc];
+        FAssertLF(memoryBuffer);
+        memset(memoryBuffer, 0, memoryToAlloc);
+
+        // Move the pointers to the correct position
+        blockHeader         = reinterpret_cast<BlockHeader*>(memoryBuffer);
+        blockIndexesTable   = reinterpret_cast<int*>(memoryBuffer+sizeof(BlockHeader));
+        blockCells          = reinterpret_cast<CellClass*>(memoryBuffer+sizeof(BlockHeader)+(blockIndexesTableSize*sizeof(int)));
+
+        // Init header
+        blockHeader->startingIndex = inStartingIndex;
+        blockHeader->endingIndex   = inEndingIndex;
+        blockHeader->numberOfCellsInBlock  = inNumberOfCells;
+        blockHeader->blockIndexesTableSize = blockIndexesTableSize;
+
+        // Set all index to not used
+        for(int idxCellPtr = 0 ; idxCellPtr < blockIndexesTableSize ; ++idxCellPtr){
+            blockIndexesTable[idxCellPtr] = CellIsEmptyFlag;
+        }
+    }
+
+    /** Call the destructor of cells and dealloc block memory */
+    ~FGroupOfCells(){
+        for(int idxCellPtr = 0 ; idxCellPtr < blockHeader->blockIndexesTableSize ; ++idxCellPtr){
+            if(blockIndexesTable[idxCellPtr] != CellIsEmptyFlag){
+                (&blockCells[blockIndexesTable[idxCellPtr]])->~CellClass();
+            }
+        }
+        delete[] memoryBuffer;
+    }
+
+    /** The index of the fist cell (set from the constructor) */
+    MortonIndex getStartingIndex() const {
+        return blockHeader->startingIndex;
+    }
+
+    /** The index of the last cell + 1 (set from the constructor) */
+    MortonIndex getEndingIndex() const {
+        return blockHeader->endingIndex;
+    }
+
+    /** The number of cell (set from the constructor) */
+    int getNumberOfCellsInBlock() const {
+        return blockHeader->numberOfCellsInBlock;
+    }
+
+    /** The size of the interval endingIndex-startingIndex (set from the constructor) */
+    int getSizeOfInterval() const {
+        return blockHeader->blockIndexesTableSize;
+    }
+
+    /** Return true if inIndex should be located in the current block */
+    bool isInside(const MortonIndex inIndex) const{
+        return blockHeader->startingIndex <= inIndex && inIndex < blockHeader->endingIndex;
+    }
+
+    /** Return true if inIndex is located in the current block and is not empty */
+    bool exists(const MortonIndex inIndex) const {
+        return isInside(inIndex) && (blockIndexesTable[inIndex-blockHeader->startingIndex] != CellIsEmptyFlag);
+    }
+
+    /** Return the address of the cell if it exists (or NULL) */
+    CellClass* getCell(const MortonIndex inIndex){
+        if( exists(inIndex) ) return &blockCells[blockIndexesTable[inIndex-blockHeader->startingIndex]];
+        else return 0;
+    }
+
+    /** Return the address of the cell if it exists (or NULL) */
+    const CellClass* getCell(const MortonIndex inIndex) const {
+        if( exists(inIndex) ) return &blockCells[blockIndexesTable[inIndex-blockHeader->startingIndex]];
+        else return 0;
+    }
+
+    /** Allocate a new cell by calling its constructor */
+    template<typename... CellConstructorParams>
+    void newCell(const MortonIndex inIndex, const int id, CellConstructorParams... args){
+        FAssertLF(isInside(inIndex));
+        FAssertLF(!exists(inIndex));
+        FAssertLF(id < blockHeader->blockIndexesTableSize);
+        new((void*)&blockCells[id]) CellClass(args...);
+        blockIndexesTable[inIndex-blockHeader->startingIndex] = id;
+    }
+
+    /** Iterate on each allocated cells */
+    template<typename... FunctionParams>
+    void forEachCell(std::function<void(CellClass*)> function, FunctionParams... args){
+        for(int idxCellPtr = 0 ; idxCellPtr < blockHeader->blockIndexesTableSize ; ++idxCellPtr){
+            if(blockIndexesTable[idxCellPtr] != CellIsEmptyFlag){
+                function(&blockCells[blockIndexesTable[idxCellPtr]], args...);
+            }
+        }
+    }
+};
+
+
+#endif // FGROUPOFCELLS_HPP

Src/GroupTree/FGroupOfParticles.hpp

+#ifndef FGROUPOFPARTICLES_HPP
+#define FGROUPOFPARTICLES_HPP
+
+
+#include "../Utils/FAssert.hpp"
+#include "../Containers/FTreeCoordinate.hpp"
+
+#include <list>
+#include <functional>
+
+
+/**
+* @brief The FGroupOfParticles class manages the leaves in block allocation.
+*/
+template <unsigned NbAttributesPerParticle, class AttributeClass = FReal>
+class FGroupOfParticles {
+    /** One header is allocated at the beginning of each block */
+    struct BlockHeader{
+        MortonIndex startingIndex;
+        MortonIndex endingIndex;
+        int numberOfLeavesInBlock;
+        int blockIndexesTableSize;
+    };
+
+    /** Information about a leaf */
+    struct LeafHeader {
+        int nbParticles;
+        size_t offSet;
+    };
+
+protected:
+    //< Pointer to a block memory
+    unsigned char* memoryBuffer;
+
+    //< Pointer to the header inside the block memory
+    BlockHeader*    blockHeader;
+    //< Pointer to the indexes table inside the block memory
+    int*            blockIndexesTable;
+    //< Pointer to leaves information
+    LeafHeader*     leafHeader;
+    //< The total number of particles in the group
+    const int nbParticlesInGroup;
+
+    //< Pointers to particle position x, y, z
+    FReal* particlePosition[3];
+    //< Bytes difference/offset between position
+    size_t positionOffset;
+
+    //< Pointers to the particles data inside the block memory
+    AttributeClass*      particleAttributes[NbAttributesPerParticle];
+    //< Bytes difference/offset between attributes
+    size_t attributeOffset;
+
+    //< This value is for not used leaves
+    static const MortonIndex LeafIsEmptyFlag = -1;
+
+public:
+    /**
+ * @brief FGroupOfParticles
+ * @param inStartingIndex first leaf morton index
+ * @param inEndingIndex last leaf morton index + 1
+ * @param inNumberOfLeaves total number of leaves in the interval (should be <= inEndingIndex-inEndingIndex)
+ */
+    FGroupOfParticles(const MortonIndex inStartingIndex, const MortonIndex inEndingIndex, const int inNumberOfLeaves, const int inNbParticles)
+        : memoryBuffer(0), blockHeader(0), blockIndexesTable(0), leafHeader(0), nbParticlesInGroup(inNbParticles),
+          positionOffset(0), attributeOffset(0) {
+        memset(particlePosition, 0, sizeof(particlePosition));
+        memset(particleAttributes, 0, sizeof(particleAttributes));
+
+        // Find the number of leaf to allocate in the blocks
+        const int blockIndexesTableSize = int(inEndingIndex-inStartingIndex);
+        FAssertLF(inNumberOfLeaves <= blockIndexesTableSize);
+        // Total number of bytes in the block
+        const size_t sizeOfOneParticle = (3*sizeof(FReal) + NbAttributesPerParticle*sizeof(AttributeClass));
+        const size_t memoryToAlloc = sizeof(BlockHeader) + (blockIndexesTableSize*sizeof(int)) + (inNumberOfLeaves*sizeof(LeafHeader))
+                + inNbParticles*sizeOfOneParticle;
+
+        // Allocate
+        memoryBuffer = new unsigned char[memoryToAlloc];
+        FAssertLF(memoryBuffer);
+        memset(memoryBuffer, 0, memoryToAlloc);
+
+        // Move the pointers to the correct position
+        blockHeader         = reinterpret_cast<BlockHeader*>(memoryBuffer);
+        blockIndexesTable   = reinterpret_cast<int*>(memoryBuffer+sizeof(BlockHeader));
+        leafHeader          = reinterpret_cast<LeafHeader*>(memoryBuffer+sizeof(BlockHeader)+(blockIndexesTableSize*sizeof(int)));
+
+        // Init header
+        blockHeader->startingIndex = inStartingIndex;
+        blockHeader->endingIndex   = inEndingIndex;
+        blockHeader->numberOfLeavesInBlock  = inNumberOfLeaves;
+        blockHeader->blockIndexesTableSize = blockIndexesTableSize;
+
+        // Init particle pointers
+        positionOffset = (sizeof(FReal) * inNbParticles);
+        particlePosition[0] = reinterpret_cast<FReal*>(leafHeader + inNumberOfLeaves);
+        particlePosition[1] = (particlePosition[0] + inNbParticles);
+        particlePosition[2] = (particlePosition[1] + inNbParticles);
+
+        // Redirect pointer to data
+        attributeOffset = (sizeof(AttributeClass) * inNbParticles);
+        unsigned char* previousPointer = reinterpret_cast<unsigned char*>(particlePosition[2] + inNbParticles);
+        for(unsigned idxAttribute = 0 ; idxAttribute < NbAttributesPerParticle ; ++idxAttribute){
+            particleAttributes[idxAttribute] = reinterpret_cast<AttributeClass*>(previousPointer);
+            previousPointer += sizeof(AttributeClass)*inNbParticles;
+        }
+
+        // Set all index to not used
+        for(int idxLeafPtr = 0 ; idxLeafPtr < blockIndexesTableSize ; ++idxLeafPtr){
+            blockIndexesTable[idxLeafPtr] = LeafIsEmptyFlag;
+        }
+    }
+
+    /** Call the destructor of leaves and dealloc block memory */
+    ~FGroupOfParticles(){
+        delete[] memoryBuffer;
+    }
+
+    /** The index of the fist leaf (set from the constructor) */
+    MortonIndex getStartingIndex() const {
+        return blockHeader->startingIndex;
+    }
+
+    /** The index of the last leaf + 1 (set from the constructor) */
+    MortonIndex getEndingIndex() const {
+        return blockHeader->endingIndex;
+    }
+
+    /** The number of leaf (set from the constructor) */
+    int getNumberOfLeavesInBlock() const {
+        return blockHeader->numberOfLeavesInBlock;
+    }
+
+    /** Get the total number of particles in the group */
+    int getNbParticlesInGroup() const {
+        return nbParticlesInGroup;
+    }
+
+    /** The size of the interval endingIndex-startingIndex (set from the constructor) */
+    int getSizeOfInterval() const {
+        return blockHeader->blockIndexesTableSize;
+    }
+
+    /** Return true if inIndex should be located in the current block */
+    bool isInside(const MortonIndex inIndex) const{
+        return blockHeader->startingIndex <= inIndex && inIndex < blockHeader->endingIndex;
+    }
+
+    /** Return true if inIndex is located in the current block and is not empty */
+    bool exists(const MortonIndex inIndex) const {
+        return isInside(inIndex) && (blockIndexesTable[inIndex-blockHeader->startingIndex] != LeafIsEmptyFlag);
+    }
+
+    /** Allocate a new leaf by calling its constructor */
+    template<class ParticleClassContainer>
+    void newLeaf(const MortonIndex inIndex, const int id, const int nbParticles, const size_t offsetInGroup,
+                 const ParticleClassContainer* particles, const int offsetInSrcContainer){
+        FAssertLF(isInside(inIndex));
+        FAssertLF(!exists(inIndex));
+        FAssertLF(id < blockHeader->blockIndexesTableSize);
+        blockIndexesTable[inIndex-blockHeader->startingIndex] = id;
+        leafHeader[id].nbParticles = nbParticles;
+        leafHeader[id].offSet = offsetInGroup;
+
+        // Copy position
+        memcpy(particlePosition[0] + offsetInGroup, particles->getPositions()[0] + offsetInSrcContainer, nbParticles*sizeof(FReal));
+        memcpy(particlePosition[1] + offsetInGroup, particles->getPositions()[1] + offsetInSrcContainer, nbParticles*sizeof(FReal));
+        memcpy(particlePosition[2] + offsetInGroup, particles->getPositions()[2] + offsetInSrcContainer, nbParticles*sizeof(FReal));
+
+        // Copy data
+        for(unsigned idxAttribute = 0 ; idxAttribute < NbAttributesPerParticle ; ++idxAttribute){
+            memcpy(particleAttributes[idxAttribute] + offsetInGroup, particles->getAttribute(idxAttribute) + offsetInSrcContainer, nbParticles*sizeof(AttributeClass));
+        }
+    }
+
+    /** Iterate on each allocated leaves */
+    template<class ParticlesAttachedClass>
+    void forEachLeaf(std::function<void(ParticlesAttachedClass*)> function){
+        for(int idxLeafPtr = 0 ; idxLeafPtr < blockHeader->blockIndexesTableSize ; ++idxLeafPtr){
+            if(blockIndexesTable[idxLeafPtr] != LeafIsEmptyFlag){
+                const int id = blockIndexesTable[idxLeafPtr];
+                ParticlesAttachedClass leaf(leafHeader[id].nbParticles,
+                                            particlePosition[0] + leafHeader[id].offSet,
+                        positionOffset,
+                        particleAttributes[0] + leafHeader[id].offSet,
+                        attributeOffset);
+                function(&leaf);
+            }
+        }
+    }
+
+
+    /** Return the address of the leaf if it exists (or NULL) */
+    template<class ParticlesAttachedClass>
+    ParticlesAttachedClass* getLeaf(const MortonIndex leafIndex){
+        if(blockIndexesTable[leafIndex - blockHeader->startingIndex] != LeafIsEmptyFlag){
+            const int id = blockIndexesTable[leafIndex - blockHeader->startingIndex];
+            return ParticlesAttachedClass(leafHeader[id].nbParticles,
+                                          particlePosition[0] + leafHeader[id].offSet,
+                    positionOffset,
+                    particleAttributes[0] + leafHeader[id].offSet,
+                    attributeOffset);
+        }
+        return NULL;
+    }
+};
+
+
+
+#endif // FGROUPOFPARTICLES_HPP

Src/GroupTree/FGroupTree.hpp

+#ifndef FGROUPTREE_HPP
+#define FGROUPTREE_HPP
+
+#include "../Utils/FAssert.hpp"
+#include "../Containers/FTreeCoordinate.hpp"
+
+#include "FGroupOfCells.hpp"
+#include "FGroupOfParticles.hpp"
+
+#include <list>
+#include <functional>
+
+template <class CellClass, unsigned NbAttributesPerParticle, class AttributeClass = FReal>
+class FGroupTree {
+    //< This value is for not used cells
+    static const int CellIsEmptyFlag = -1;
+
+protected:
+    //< 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::list<FGroupOfCells<CellClass>*>* cellBlocksPerLevel;
+    //< all the blocks of leaves
+    std::list<FGroupOfParticles<NbAttributesPerParticle,AttributeClass>*> particleBlocks;
+
+    /** @brief This private method take an array of Morton index to
+   * create the block and the cells, using the constructor of
+   * FGroupOfCells !!
+   */
+    FGroupOfCells<CellClass> * createBlockFromArray(MortonIndex head[]){
+        //Store the start and end
+        MortonIndex start = head[0];
+        MortonIndex end = start;
+        int count = 0;
+        // Find the number of cell to allocate in the blocks
+        for(int idxHead = 0 ; idxHead<nbElementsPerBlock ; idxHead++){
+            if (head[idxHead] != CellIsEmptyFlag){
+                count++;
+                end = head[idxHead];
+            }
+            // else{
+            // 	break;
+            // }
+        }
+        //allocation of memory
+        FGroupOfCells<CellClass> * newBlock = new FGroupOfCells<CellClass>(start,end+1,count);
+        //allocation of cells
+        for(int idx=0 ; idx<count ; idx++){
+            newBlock->newCell(head[idx], idx);
+            (newBlock->getCell(head[idx]))->setMortonIndex(head[idx]);
+            //(this->getCell(head[idx]))->setCoordinate();
+        }
+        return newBlock;
+    }
+
+public:
+    /** This constructor create a blocked octree from a usual octree
+   * 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.
+   */
+    template<class OctreeClass>
+    FGroupTree(const int inTreeHeight, const int inNbElementsPerBlock, OctreeClass*const inOctreeSrc)
+        : treeHeight(inTreeHeight), nbElementsPerBlock(inNbElementsPerBlock), cellBlocksPerLevel(0){
+        cellBlocksPerLevel = new std::list<FGroupOfCells<CellClass>*>[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;
+                int nbParticlesInGroup = octreeIterator.getCurrentLeaf()->getSrc()->getNbParticles();
+                while(sizeOfBlock < nbElementsPerBlock && octreeIterator.moveRight()){
+                    sizeOfBlock += 1;
+                    nbParticlesInGroup += octreeIterator.getCurrentLeaf()->getSrc()->getNbParticles();
+                }
+
+                // Create a block with the apropriate parameters
+                FGroupOfCells<CellClass>*const newBlock = new FGroupOfCells<CellClass>(blockIteratorInOctree.getCurrentGlobalIndex(),
+                                                                 octreeIterator.getCurrentGlobalIndex()+1,
+                                                                 sizeOfBlock);
+                FGroupOfParticles<NbAttributesPerParticle, AttributeClass>*const newParticleBlock = new FGroupOfParticles<NbAttributesPerParticle, AttributeClass>(blockIteratorInOctree.getCurrentGlobalIndex(),
+                                                                 octreeIterator.getCurrentGlobalIndex()+1,
+                                                                 sizeOfBlock, nbParticlesInGroup);
+
+                // Initialize each cell of the block
+                int cellIdInBlock = 0;
+                int nbParticlesBeforeLeaf = 0;
+                while(cellIdInBlock != sizeOfBlock){
+                    // Add cell
+                    const CellClass*const oldNode = blockIteratorInOctree.getCurrentCell();
+                    newBlock->newCell(oldNode->getMortonIndex(), cellIdInBlock);
+
+                    CellClass* newNode = newBlock->getCell(oldNode->getMortonIndex());
+                    newNode->setMortonIndex(oldNode->getMortonIndex());
+                    newNode->setCoordinate(oldNode->getCoordinate());
+
+                    // Add leaf
+                    newParticleBlock->newLeaf(oldNode->getMortonIndex(), cellIdInBlock,
+                                              blockIteratorInOctree.getCurrentLeaf()->getSrc()->getNbParticles(),
+                                              nbParticlesBeforeLeaf, blockIteratorInOctree.getCurrentLeaf()->getSrc(), 0);
+
+                    nbParticlesBeforeLeaf += blockIteratorInOctree.getCurrentLeaf()->getSrc()->getNbParticles();
+                    cellIdInBlock += 1;
+                    blockIteratorInOctree.moveRight();
+                }
+
+                // Keep the block
+                cellBlocksPerLevel[idxLevel].push_back(newBlock);
+                particleBlocks.push_back(newParticleBlock);
+
+                // If we can move right then add another block
+            } while(octreeIterator.moveRight());
+
+            avoidGotoLeft.moveUp();
+            octreeIterator = avoidGotoLeft;
+        }
+
+        // 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;