diff --git a/Src/GroupTree/Core/FGroupTree.hpp b/Src/GroupTree/Core/FGroupTree.hpp index b1f4ba42c515dd7db381c9e8f798bb74db8d448d..8c8c4269e055d8616b509378e6ebe5d039270725 100644 --- a/Src/GroupTree/Core/FGroupTree.hpp +++ b/Src/GroupTree/Core/FGroupTree.hpp @@ -169,184 +169,6 @@ public: } } - - /** - * 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. - */ - template<class ParticleContainer> - FGroupTree(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, - const int inNbElementsPerBlock, ParticleContainer* inParticlesContainer, const bool particlesAreSorted = false): - 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; - }); - } - - // 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){ - CellGroupConstIterator iterChildCells = cellBlocksPerLevel[idxLevel+1].begin(); - const CellGroupConstIterator iterChildEndCells = cellBlocksPerLevel[idxLevel+1].end(); - - MortonIndex currentCellIndex = (*iterChildCells)->getStartingIndex(); - 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(); - } - } - } - - // 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; - } - } - } - 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. @@ -358,7 +180,7 @@ public: template<class ParticleContainer> FGroupTree(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const int inNbElementsPerBlock, ParticleContainer* inParticlesContainer, - const bool particlesAreSorted, MortonIndex inLeftLimite): + const bool particlesAreSorted = false, MortonIndex inLeftLimite = -1): treeHeight(inTreeHeight),nbElementsPerBlock(inNbElementsPerBlock),cellBlocksPerLevel(nullptr), boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth), boxWidthAtLeafLevel(inBoxWidth/FReal(1<<(inTreeHeight-1))){