Add the possibility to create one parent block for each subblock

Src/GroupTree/Core/FGroupTree.hpp

@@ -362,6 +362,266 @@ public:
         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 inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter,
+               const int inNbElementsPerBlock, ParticleContainer* inParticlesContainer,
+               const bool particlesAreSorted, const bool oneParent,
+               const FReal inCoverRatio = 0.0, MortonIndex inLeftLimite = -1):
+            treeHeight(inTreeHeight),nbElementsPerBlock(inNbElementsPerBlock),cellBlocksPerLevel(nullptr),
+            boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth),
+            boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(inTreeHeight-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{
+                int originalIndex;
+                MortonIndex mindex;
+            };
+            // Convert position to morton index
+            const int 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(int 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(treeHeight-1);
+                    particlesToSort[idxPart].mindex = particleIndex;
+                    particlesToSort[idxPart].originalIndex = idxPart;
+                }
+            }
+
+            // Sort if needed
+            if(particlesAreSorted == false){
+                FQuickSort<ParticleSortingStruct, int>::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){
+                    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);
+
+                    CompositeCellClass newNode = newBlock->getCompleteCell(currentBlockIndexes[cellIdInBlock]);
+                    newNode.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
+                    FTreeCoordinate coord;
+                    coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock], idxLevel);
+                    newNode.setCoordinate(coord);
+
+                    // Add leaf
+                    nbParticlesOffsetBeforeLeaf = newParticleBlock->newLeaf(currentBlockIndexes[cellIdInBlock], cellIdInBlock,
+                                              nbParticlesPerLeaf[cellIdInBlock], nbParticlesOffsetBeforeLeaf);
+
+                    BasicAttachedClass attachedLeaf = newParticleBlock->template getLeaf<BasicAttachedClass>(currentBlockIndexes[cellIdInBlock]);
+                    // Copy each particle from the original position
+                    for(int idxPart = 0 ; idxPart < nbParticlesPerLeaf[cellIdInBlock] ; ++idxPart){
+                        attachedLeaf.setParticle(idxPart, particlesToSort[idxPart + offsetParticles].originalIndex, inParticlesContainer);
+                    }
+                    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;
+
+            if(oneParent == false){
+                // We need to proceed each group in sub level
+                while(iterChildCells != iterChildEndCells){
+                    // Count until end of sub group is reached or we have enough cells
+                    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();
+                            }
+                        }
+                    }
+
+                    // If group is full
+                    if(sizeOfBlock == nbElementsPerBlock || (sizeOfBlock && iterChildCells == iterChildEndCells)){
+                        // Create a group
+                        CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
+                                                                         currentBlockIndexes[sizeOfBlock-1]+1,
+                                                                         sizeOfBlock);
+                        // Init cells
+                        for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
+                            newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);
+
+                            CompositeCellClass newNode = newBlock->getCompleteCell(currentBlockIndexes[cellIdInBlock]);
+                            newNode.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
+                            FTreeCoordinate coord;
+                            coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock], idxLevel);
+                            newNode.setCoordinate(coord);
+                        }
+
+                        // Keep the block
+                        cellBlocksPerLevel[idxLevel].push_back(newBlock);
+
+                        sizeOfBlock = 0;
+                    }
+                }
+            }
+            else{
+                // We need to proceed each group in sub level
+                while(iterChildCells != iterChildEndCells){
+                    // We want one parent group per child group so we will stop the parent group
+                    // when we arrive to the same parent as lastChildIndex (which is lastChildIndex>>3)
+                    const MortonIndex lastChildIndex = ((*iterChildCells)->getEndingIndex()-1);
+                    // Count until end of sub group is reached or we passe the requested parent
+                    while( iterChildCells != iterChildEndCells
+                           && (currentCellIndex>>3) <= (lastChildIndex>>3) ){
+                        // Proceed until the requested parent
+                        while(currentCellIndex != (*iterChildCells)->getEndingIndex()
+                              && (currentCellIndex>>3) <= (lastChildIndex>>3) ){
+                            if((*iterChildCells)->exists(currentCellIndex)){
+                                currentBlockIndexes[sizeOfBlock] = (currentCellIndex>>3);
+                                sizeOfBlock += 1;
+                                currentCellIndex = (((currentCellIndex>>3)+1)<<3);
+                            }
+                            else{
+                                currentCellIndex += 1;
+                            }
+                        }
+                        // If we are at the end of the sub group, move to next (otherwise we have consume a part of it)
+                        if((*iterChildCells)->getEndingIndex() <= currentCellIndex){
+                            ++iterChildCells;
+                            // Update morton index
+                            if(iterChildCells != iterChildEndCells && currentCellIndex < (*iterChildCells)->getStartingIndex()){
+                                currentCellIndex = (*iterChildCells)->getStartingIndex();
+                            }
+                        }
+                    }
+
+                    // If group is full
+                    if(sizeOfBlock){
+                        // Create a group
+                        CellGroupClass*const newBlock = new CellGroupClass(currentBlockIndexes[0],
+                                                                         currentBlockIndexes[sizeOfBlock-1]+1,
+                                                                         sizeOfBlock);
+                        // Init cells
+                        for(int cellIdInBlock = 0; cellIdInBlock != sizeOfBlock ; ++cellIdInBlock){
+                            newBlock->newCell(currentBlockIndexes[cellIdInBlock], cellIdInBlock);
+
+                            CompositeCellClass newNode = newBlock->getCompleteCell(currentBlockIndexes[cellIdInBlock]);
+                            newNode.setMortonIndex(currentBlockIndexes[cellIdInBlock]);
+                            FTreeCoordinate coord;
+                            coord.setPositionFromMorton(currentBlockIndexes[cellIdInBlock], idxLevel);
+                            newNode.setCoordinate(coord);
+                        }
+
+                        // Keep the block
+                        cellBlocksPerLevel[idxLevel].push_back(newBlock);
+
+                        sizeOfBlock = 0;
+                    }
+                    else{
+                        assert(iterChildCells == iterChildEndCells);
+                    }
+                }
+            }
+        }
+        delete[] currentBlockIndexes;
+    }
+
+
     /** This function dealloc the tree by deleting each block */
     ~FGroupTree(){
         for(int idxLevel = 0 ; idxLevel < treeHeight ; ++idxLevel){

Tests/noDist/testBlockedAlgorithm.cpp

@@ -122,6 +122,7 @@ int main(int argc, char* argv[]){
     // Put the data into the tree
     //GroupOctreeClass groupedTree(NbLevels, groupSize, &tree);
     GroupOctreeClass groupedTree(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(), groupSize, &allParticles);
+    //GroupOctreeClass groupedTree(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(), groupSize, &allParticles, false, true);
     groupedTree.printInfoBlocks();
 
     // Check tree structure at leaf level