Update the way of computing mpi

Src/Components/FAbstractSendable.hpp

@@ -43,6 +43,8 @@ protected:
         static_assert(sizeof(BufferWriterClass) == 0 , "Your class should implement deserializeUp");
     }
 
+    virtual int getSavedSizeUp() = 0;
+
     ///////////////////////////////////////////////
     // For Downward pass
     ///////////////////////////////////////////////
@@ -57,6 +59,8 @@ protected:
     void deserializeDown(BufferReaderClass&){
         static_assert(sizeof(BufferWriterClass) == 0 , "Your class should implement deserializeDown");
     }
+
+    virtual int getSavedSizeDown() = 0;
 };
 
 

Src/Components/FAbstractSerializable.hpp

@@ -45,14 +45,7 @@ protected:
         static_assert(sizeof(BufferReaderClass) == 0 , "Your class should implement restore");
     }
 
-    /**
-   * @brief GetSize returns the size of the current block (in order to prepare the buffer)
-   * @return
-   */
-    static constexpr int GetSize(){
-        return 0;
-    }
-
+    virtual int getSavedSize() const = 0;
 };
 
 #endif // FABSTRACTSERIALIZABLE_HPP

Src/Components/FBasicCell.hpp

@@ -53,6 +53,10 @@ public:
         FExtendCoordinate::restore(buffer);
     }
 
+    int getSavedSize() const {
+        return FExtendMortonIndex::getSavedSize() + FExtendMortonIndex::getSavedSize();
+    }
+
     /** Do nothing */
     void resetToInitialState(){
     }

Src/Components/FTestCell.hpp

@@ -76,10 +76,9 @@ public:
 		buffer >> dataDown >> dataUp;
 	}
 
-	//   	  static constexpr int GetSize(){
-	static constexpr     std::size_t GetSize(){
-		return sizeof(long long int)*2;
-	}
+    int getSavedSize() const {
+        return int(sizeof(long long int))*2;
+    }
 
 	/////////////////////////////////////////////////
 
@@ -105,6 +104,13 @@ public:
 		buffer >> this->dataDown;
 	}
 
+    int getSavedSizeDown() {
+        return int(sizeof(long long int));
+    }
+
+    int getSavedSizeUp() {
+        return int(sizeof(long long int));
+    }
 };
 
 

Src/Containers/FTreeCoordinate.hpp

@@ -258,6 +258,11 @@ public:
         buffer >> data[0] >> data[1] >> data[2];
     }
 
+    /** To know the size when we save it */
+    int getSavedSize() const {
+        return int(sizeof(data[0]) + sizeof(data[1]) + sizeof(data[2]));
+    }
+
 
     static std::string MortonToBinary(MortonIndex index, int level){
         std::string str;

Src/Core/FFmmAlgorithmThreadProc.hpp

@@ -62,9 +62,6 @@
  */
 template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
 class FFmmAlgorithmThreadProc : public FAbstractAlgorithm {
-    // Can be deleted
-    // const static int MaxSizePerCell = CellClass::GetSize();
-
     OctreeClass* const tree;                 //< The octree to work on
     KernelClass** kernels;                   //< The kernels
 
@@ -533,7 +530,6 @@ private:
 
 
     void transferPass(){
-        const int MaxSizePerCell = CellClass::GetSize();
         FLOG( FLog::Controller.write("\tStart Downward Pass (M2L)\n").write(FLog::Flush); );
         FLOG(FTic counterTime);
         FLOG(FTic computationCounter);
@@ -549,16 +545,16 @@ private:
         // pointer to send
         FVector<typename OctreeClass::Iterator> toSend[nbProcess * OctreeHeight];
         // index
-        int*const indexToSend = new int[nbProcess * OctreeHeight];
-        memset(indexToSend, 0, sizeof(int) * nbProcess * OctreeHeight);
+        long long int*const indexToSend = new long long int[nbProcess * OctreeHeight];
+        memset(indexToSend, 0, sizeof(long long int) * nbProcess * OctreeHeight);
         // To know which one has need someone
         FBoolArray** const leafsNeedOther = new FBoolArray*[OctreeHeight];
         memset(leafsNeedOther, 0, sizeof(FBoolArray*) * OctreeHeight);
 
         // All process say to each others
         // what the will send to who
-        int*const globalReceiveMap = new int[nbProcess * nbProcess * OctreeHeight];
-        memset(globalReceiveMap, 0, sizeof(int) * nbProcess * nbProcess * OctreeHeight);
+        long long int*const globalReceiveMap = new long long  int[nbProcess * nbProcess * OctreeHeight];
+        memset(globalReceiveMap, 0, sizeof(long long  int) * nbProcess * nbProcess * OctreeHeight);
 
         FMpiBufferWriter**const sendBuffer = new FMpiBufferWriter*[nbProcess * OctreeHeight];
         memset(sendBuffer, 0, sizeof(FMpiBufferWriter*) * nbProcess * OctreeHeight);
@@ -566,9 +562,9 @@ private:
         FMpiBufferReader**const recvBuffer = new FMpiBufferReader*[nbProcess * OctreeHeight];
         memset(recvBuffer, 0, sizeof(FMpiBufferReader*) * nbProcess * OctreeHeight);
 
-#pragma omp parallel
+        #pragma omp parallel
         {
-#pragma omp master
+            #pragma omp master
             {
                 {
                     FLOG(prepareCounter.tic());
@@ -636,7 +632,11 @@ private:
                                         needOther = true;
 
                                         toSend[idxLevel * nbProcess + procToReceive].push(iterArrayLocal[idxCell]);
-                                        ++indexToSend[idxLevel * nbProcess + procToReceive];
+                                        if(indexToSend[idxLevel * nbProcess + procToReceive] == 0){
+                                            indexToSend[idxLevel * nbProcess + procToReceive] = sizeof(int);
+                                        }
+                                        indexToSend[idxLevel * nbProcess + procToReceive] += iterArrayLocal[idxCell].getCurrentCell()->getSavedSizeUp();
+                                        indexToSend[idxLevel * nbProcess + procToReceive] += sizeof(MortonIndex);
                                     }
                                 }
                             }
@@ -655,7 +655,7 @@ private:
                 //////////////////////////////////////////////////////////////////
 
                 FLOG(gatherCounter.tic());
-                FMpi::MpiAssert( MPI_Allgather( indexToSend, nbProcess * OctreeHeight, MPI_INT, globalReceiveMap, nbProcess * OctreeHeight, MPI_INT, comm.getComm()),  __LINE__ );
+                FMpi::MpiAssert( MPI_Allgather( indexToSend, nbProcess * OctreeHeight, MPI_LONG_LONG_INT, globalReceiveMap, nbProcess * OctreeHeight, MPI_LONG_LONG_INT, comm.getComm()),  __LINE__ );
                 FLOG(gatherCounter.tac());
 
                 //////////////////////////////////////////////////////////////////
@@ -669,28 +669,30 @@ private:
                 MPI_Status*const status = new MPI_Status[2 * nbProcess * OctreeHeight];
                 int iterRequest = 0;
 
-                const int SizeOfCellToSend = sizeof(MortonIndex) + sizeof(int) + MaxSizePerCell;
-
                 for(int idxLevel = 2 ; idxLevel < OctreeHeight ; ++idxLevel ){
                     for(int idxProc = 0 ; idxProc < nbProcess ; ++idxProc){
-                        const int toSendAtProcAtLevel = indexToSend[idxLevel * nbProcess + idxProc];
+                        const long long int toSendAtProcAtLevel = indexToSend[idxLevel * nbProcess + idxProc];
                         if(toSendAtProcAtLevel != 0){
-                            sendBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferWriter(comm.getComm(),toSendAtProcAtLevel * SizeOfCellToSend);
+                            sendBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferWriter(comm.getComm(),toSendAtProcAtLevel);
 
-                            for(int idxLeaf = 0 ; idxLeaf < toSendAtProcAtLevel; ++idxLeaf){
+                            sendBuffer[idxLevel * nbProcess + idxProc]->write(int(toSend[idxLevel * nbProcess + idxProc].getSize()));
+
+                            for(int idxLeaf = 0 ; idxLeaf < toSend[idxLevel * nbProcess + idxProc].getSize(); ++idxLeaf){
                                 const MortonIndex cellIndex = toSend[idxLevel * nbProcess + idxProc][idxLeaf].getCurrentGlobalIndex();
                                 sendBuffer[idxLevel * nbProcess + idxProc]->write(cellIndex);
                                 toSend[idxLevel * nbProcess + idxProc][idxLeaf].getCurrentCell()->serializeUp(*sendBuffer[idxLevel * nbProcess + idxProc]);
                             }
 
+                            FAssertLF(sendBuffer[idxLevel * nbProcess + idxProc]->getSize() == toSendAtProcAtLevel);
+
                             FMpi::MpiAssert( MPI_Isend( sendBuffer[idxLevel * nbProcess + idxProc]->data(),
                                              sendBuffer[idxLevel * nbProcess + idxProc]->getSize(),MPI_PACKED, idxProc,
                                     FMpi::TagLast + idxLevel, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
                         }
 
-                        const int toReceiveFromProcAtLevel = globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess];
+                        const long long int toReceiveFromProcAtLevel = globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess];
                         if(toReceiveFromProcAtLevel){
-                            recvBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferReader(comm.getComm(),toReceiveFromProcAtLevel * SizeOfCellToSend);
+                            recvBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferReader(comm.getComm(),toReceiveFromProcAtLevel);
 
                             FMpi::MpiAssert( MPI_Irecv(recvBuffer[idxLevel * nbProcess + idxProc]->data(),
                                              recvBuffer[idxLevel * nbProcess + idxProc]->getCapacity(), MPI_PACKED,idxProc,
@@ -704,7 +706,7 @@ private:
                 //////////////////////////////////////////////////////////////////
 
                 // Wait to receive every things (and send every things)
-                MPI_Waitall(iterRequest, requests, status);
+                FMpi::MpiAssert(MPI_Waitall(iterRequest, requests, status), __LINE__);
 
                 delete[] requests;
                 delete[] status;
@@ -716,7 +718,7 @@ private:
             // Do M2L
             //////////////////////////////////////////////////////////////////
 
-#pragma omp single nowait
+            #pragma omp single nowait
             {
                 typename OctreeClass::Iterator octreeIterator(tree);
                 octreeIterator.moveDown();
@@ -746,7 +748,7 @@ private:
                     {
                         const int chunckSize = FMath::Max(1, numberOfCells/(omp_get_num_threads()*omp_get_num_threads()));
                         for(int idxCell = 0 ; idxCell < numberOfCells ; idxCell += chunckSize){
-#pragma omp task default(none) shared(numberOfCells,idxLevel) firstprivate(idxCell) //+ shared(chunckSize)
+                            #pragma omp task default(none) shared(numberOfCells,idxLevel) firstprivate(idxCell) //+ shared(chunckSize)
                             {
                                 KernelClass * const myThreadkernels = kernels[omp_get_thread_num()];
                                 const CellClass* neighbors[343];
@@ -760,15 +762,15 @@ private:
                         }
                     }//End of task spawning
 
-#pragma omp taskwait
+                    #pragma omp taskwait
 
                     for(int idxThread = 0 ; idxThread < omp_get_num_threads() ; ++idxThread){
-#pragma omp task  default(none) firstprivate(idxThread) shared(idxLevel)
+                        #pragma omp task  default(none) firstprivate(idxThread) shared(idxLevel)
                         {
                             kernels[idxThread]->finishedLevelM2L(idxLevel);
                         }
                     }
-#pragma omp taskwait
+                    #pragma omp taskwait
 
                     FLOG(computationCounter.tac());
                 }
@@ -793,20 +795,24 @@ private:
                 // put the received data into a temporary tree
                 FLightOctree<CellClass> tempTree;
                 for(int idxProc = 0 ; idxProc < nbProcess ; ++idxProc){
-                    const int toReceiveFromProcAtLevel = globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess];
+                    if(recvBuffer[idxLevel * nbProcess + idxProc]){
+                        const int toReceiveFromProcAtLevel = recvBuffer[idxLevel * nbProcess + idxProc]->FMpiBufferReader::getValue<int>();
 
-                    for(int idxCell = 0 ; idxCell < toReceiveFromProcAtLevel ; ++idxCell){
-                        const MortonIndex cellIndex = recvBuffer[idxLevel * nbProcess + idxProc]->FMpiBufferReader::getValue<MortonIndex>();
+                        for(int idxCell = 0 ; idxCell < toReceiveFromProcAtLevel ; ++idxCell){
+                            const MortonIndex cellIndex = recvBuffer[idxLevel * nbProcess + idxProc]->FMpiBufferReader::getValue<MortonIndex>();
 
-                        CellClass* const newCell = new CellClass;
-                        newCell->setMortonIndex(cellIndex);
-                        newCell->deserializeUp(*recvBuffer[idxLevel * nbProcess + idxProc]);
+                            CellClass* const newCell = new CellClass;
+                            newCell->setMortonIndex(cellIndex);
+                            newCell->deserializeUp(*recvBuffer[idxLevel * nbProcess + idxProc]);
 
-                        tempTree.insertCell(cellIndex, idxLevel, newCell);
+                            tempTree.insertCell(cellIndex, idxLevel, newCell);
+                        }
+
+                        FAssertLF(globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess] ==
+                                recvBuffer[idxLevel * nbProcess + idxProc]->tell());
                     }
                 }
 
-
                 // take cells from our octree only if they are
                 // linked to received data
                 int numberOfCells = 0;
@@ -830,14 +836,14 @@ private:
 
                 // Compute this cells
                 FLOG(computationCounter.tic());
-#pragma omp parallel
+                #pragma omp parallel
                 {
                     KernelClass * const myThreadkernels = kernels[omp_get_thread_num()];
                     MortonIndex neighborsIndex[189];
                     int neighborsPosition[189];
                     const CellClass* neighbors[343];
 
-#pragma omp for schedule(static) nowait
+                    #pragma omp for schedule(static) nowait
                     for(int idxCell = 0 ; idxCell < numberOfCells ; ++idxCell){
                         // compute indexes
                         memset(neighbors, 0, 343 * sizeof(CellClass*));
@@ -899,7 +905,6 @@ private:
     //////////////////////////////////////////////////////////////////
 
     void downardPass(){ // second L2L
-        const int MaxSizePerCell = CellClass::GetSize();
         FLOG( FLog::Controller.write("\tStart Downward Pass (L2L)\n").write(FLog::Flush); );
         FLOG(FTic counterTime);
         FLOG(FTic computationCounter);
@@ -919,8 +924,8 @@ private:
 
         const int heightMinusOne = OctreeHeight - 1;
 
-        FMpiBufferWriter sendBuffer(comm.getComm(),MaxSizePerCell);
-        FMpiBufferReader recvBuffer(comm.getComm(),MaxSizePerCell);
+        FMpiBufferWriter sendBuffer(comm.getComm());
+        FMpiBufferReader recvBuffer(comm.getComm());
 
         int righestProcToSendTo   = nbProcess - 1;
 

Src/Extensions/FExtendCellType.hpp

@@ -86,6 +86,10 @@ public:
         containsTargets = false;
         containsSources = false;
     }
+
+    int getSavedSize() const {
+        return sizeof(containsTargets) + sizeof(containsSources);
+    }
 };
 
 

Src/Extensions/FExtendCoordinate.hpp

@@ -75,6 +75,9 @@ public:
         coordinate.restore(buffer);
     }
   
+    int getSavedSize() const {
+        return coordinate.getSavedSize();
+    }
 };
 
 

Src/Extensions/FExtendMortonIndex.hpp

@@ -66,6 +66,10 @@ public:
     void restore(BufferReaderClass& buffer) {
         buffer >> mortonIndex;
     }
+
+    int getSavedSize() const {
+        return int(sizeof(mortonIndex));
+    }
 };
 
 

Src/Kernels/Chebyshev/FChebCell.hpp

@@ -116,9 +116,18 @@ public:
 		buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
 	}
 
-	static constexpr int GetSize(){
-		return int(sizeof(FReal)) * VectorSize*(NRHS+NLHS)*NVALS;
-	}
+    int getSavedSize() const {
+        return int(sizeof(FReal)) * VectorSize*(NRHS+NLHS)*NVALS + FBasicCell::getSavedSize();
+    }
+
+    int getSavedSizeUp() const {
+        return int(sizeof(FReal)) * VectorSize*(NRHS)*NVALS;
+    }
+
+    int getSavedSizeDown() const {
+        return int(sizeof(FReal)) * VectorSize*(NLHS)*NVALS;
+    }
+
 //	template <class StreamClass>
 //	const void print(StreamClass& output) const{
  template <class StreamClass>

Src/Kernels/Interpolation/FInterpCell.hpp

@@ -174,6 +174,21 @@ public:
     return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
   }
 
+  int getSavedSize() const {
+      return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal)
+              + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex)
+              + FBasicCell::getSavedSize();
+  }
+
+  int getSavedSizeUp() const {
+      return (NRHS)*NVALS*VectorSize * (int) sizeof(FReal)
+              + (NRHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
+  }
+
+  int getSavedSizeDown() const {
+      return (NLHS)*NVALS*VectorSize * (int) sizeof(FReal)
+              + (NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
+  }
 
 };
 

Src/Kernels/Rotation/FRotationCell.hpp

@@ -92,10 +92,10 @@ public:
         return local_exp;
     }
 
-  const int getArraySize() const
-  {
-    return MultipoleSize;
-  }
+    const int getArraySize() const
+    {
+        return MultipoleSize;
+    }
 
 
     /** Make it like the begining */
@@ -144,9 +144,19 @@ public:
         buffer.fillArray(multipole_exp, MultipoleSize);
         buffer.fillArray(local_exp, LocalSize);
     }
-  static constexpr int GetSize(){
-    return ((int) sizeof(FComplex)) * (MultipoleSize + LocalSize);
-  }
+
+    int getSavedSize() const {
+        return ((int) sizeof(FComplex)) * (MultipoleSize + LocalSize)
+                + FBasicCell::getSavedSize();
+    }
+
+    int getSavedSizeUp() const {
+        return ((int) sizeof(FComplex)) * (MultipoleSize);
+    }
+
+    int getSavedSizeDown() const {
+        return ((int) sizeof(FComplex)) * (LocalSize);
+    }
 };
 
 template <int P>
@@ -166,6 +176,10 @@ public:
         FRotationCell<P>::resetToInitialState();
         FExtendCellType::resetToInitialState();
     }
+
+    int getSavedSize() const {
+        return FExtendCellType::getSavedSize() + FRotationCell<P>::getSavedSize();
+    }
 };
 
 #endif // FROTATIONCELL_HPP

Src/Kernels/Spherical/FSphericalCell.hpp

@@ -151,10 +151,19 @@ public:
         buffer.fillArray(multipole_exp, PoleSize);
         buffer.fillArray(local_exp, LocalSize);
     }
-  
-  static  int GetSize(){
-    return (int) sizeof(FComplex) * (PoleSize+LocalSize);
-  }
+
+    int getSavedSize() const {
+        return ((int) sizeof(FComplex)) * (PoleSize+LocalSize)
+                + FBasicCell::getSavedSize();
+    }
+
+    int getSavedSizeUp() const {
+        return ((int) sizeof(FComplex)) * (PoleSize);
+    }
+
+    int getSavedSizeDown() const {
+        return ((int) sizeof(FComplex)) * (LocalSize);
+    }
 };
 
 int FSphericalCell::DevP(-1);
@@ -181,6 +190,10 @@ public:
         FSphericalCell::resetToInitialState();
         FExtendCellType::resetToInitialState();
     }
+
+    int getSavedSize() const {
+        return FExtendCellType::getSavedSize() + FSphericalCell::getSavedSize();
+    }
 };
 
 

Src/Kernels/Uniform/FUnifCell.hpp

@@ -165,9 +165,18 @@ public:
     buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
     buffer.fillArray(transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
   }
-  
-  static constexpr int GetSize(){
-    return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
+
+  int getSavedSize() const {
+      return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex)
+              + FBasicCell::getSavedSize();
+  }
+
+  int getSavedSizeUp() const {
+      return (NRHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
+  }
+
+  int getSavedSizeDown() const {
+      return (NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
   }
 
  template <class StreamClass>

Tests/Kernels/testSphericalGalaxyCsv.cpp

@@ -40,10 +40,11 @@
 #include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"
 
 #include "../../Src/Utils/FParameterNames.hpp"
+#include "../../Src/Arranger/FAbstractMover.hpp"
 
-class VelocityContainer : public FP2PParticleContainer<> {
-  typedef FP2PParticleContainer<> Parent;
 
+class VelocityContainer : public FP2PParticleContainer<> {
+    typedef FP2PParticleContainer<> Parent;
     FVector<FPoint> velocities;
 
 public:
@@ -84,53 +85,57 @@ public:
     }
 };
 
-struct TestParticle{
-    FPoint position;
-    FReal physicalValue;
-    FReal forces[3];
-    FReal potential;
-    FPoint velocity;
-    const FPoint& getPosition(){
-        return position;
-    }
-};
+template<class OctreeClass>
+class GalaxyMover : public FAbstractMover<OctreeClass, VelocityContainer>{
+private:
+    VelocityContainer toStoreRemovedParts;
 
-template <class ParticleClass>
-class Converter {
 public:
-    template <class ContainerClass>
-    static ParticleClass GetParticleAndRemove(ContainerClass* containers, const int idxExtract){
-        const FReal*const positionsX = containers->getPositions()[0];
-        const FReal*const positionsY = containers->getPositions()[1];
-        const FReal*const positionsZ = containers->getPositions()[2];
-        const FReal*const forcesX = containers->getForcesX();
-        const FReal*const forcesY = containers->getForcesY();
-        const FReal*const forcesZ = containers->getForcesZ();
-        const FReal*const physicalValues = containers->getPhysicalValues();
-        const FReal*const potentials = containers->getPotentials();
-        FVector<FPoint> velocites = containers->getVelocities();
-
-        TestParticle part;
-        part.position.setPosition( positionsX[idxExtract],positionsY[idxExtract],positionsZ[idxExtract]);
-        part.physicalValue = physicalValues[idxExtract];
-        part.forces[0] = forcesX[idxExtract];
-        part.forces[1] = forcesY[idxExtract];
-        part.forces[2] = forcesZ[idxExtract];
-        part.potential = potentials[idxExtract];
-        part.velocity  = velocites[idxExtract];
-
-        containers->removeParticles(&idxExtract, 1);
-
-        return part;
+    GalaxyMover() {
+    }
+
+    virtual ~GalaxyMover(){
     }
 
-    template <class OctreeClass>
-    static void Insert(OctreeClass* tree, const ParticleClass& part){
-        tree->insert(part.position , part.velocity, part.physicalValue, part.forces[0],
-                part.forces[1],part.forces[2],part.potential);
+    /** To get the position of the particle at idx idxPart in leaf lf */
+    void getParticlePosition(VelocityContainer* lf, const int idxPart, FPoint* particlePos){
+        (*particlePos) = FPoint(lf->getPositions()[0][idxPart],lf->getPositions()[1][idxPart],lf->getPositions()[2][idxPart]);
+    }
+
+    /** Remove a particle but keep it to reinsert it later*/
+    void removeFromLeafAndKeep(VelocityContainer* lf, const FPoint& particlePos, const int idxPart){
+        std::array<typename VelocityContainer::AttributesClass, VelocityContainer::NbAttributes> particleValues;
+        for(int idxAttr = 0 ; idxAttr < VelocityContainer::NbAttributes ; ++idxAttr){
+            particleValues[idxAttr] = lf->getAttribute(idxAttr)[idxPart];
+        }
+
+        toStoreRemovedParts.push(particlePos,lf->getVelocities()[idxPart],particleValues);
+
+        lf->getVelocities().removeOne(idxPart);
+        lf->removeParticles(&idxPart,1);
+    }
+
+    /** Reinsert the previously saved particles */
+    void insertAllParticles(OctreeClass* tree){
+        std::array<typename VelocityContainer::AttributesClass, VelocityContainer::NbAttributes> particleValues;
+
+        for(int idxToInsert = 0; idxToInsert<toStoreRemovedParts.getNbParticles() ; ++idxToInsert){
+            for(int idxAttr = 0 ; idxAttr < VelocityContainer::NbAttributes ; ++idxAttr){
+                particleValues[idxAttr] = toStoreRemovedParts.getAttribute(idxAttr)[idxToInsert];
+            }
+            const FPoint particlePos(toStoreRemovedParts.getPositions()[0][idxToInsert],
+                                     toStoreRemovedParts.getPositions()[1][idxToInsert],
+                                     toStoreRemovedParts.getPositions()[2][idxToInsert]);
+