Arranger works now with periodic movements of particles

Src/Arranger/FAbstractLeafInterface.hpp

@@ -21,7 +21,7 @@ class FAbstractLeafInterface{
 public:
     //    virtual FPoint getParticlePosition(ParticleClass* lf, const int idxPart) = 0;
     virtual void getParticlePosition(ParticleClass* lf, const int idxPart, FPoint& newPos) = 0;
-    virtual void removeFromLeafAndKeep(ParticleClass* lf, const int idxPart) = 0;
+    virtual void removeFromLeafAndKeep(ParticleClass* lf, const FPoint& newPos,const int idxPart) = 0;
     virtual void insertAllParticles(OctreeClass* tree) = 0;
 };
 

Src/Arranger/FArrangerPeriodic.hpp

@@ -35,7 +35,7 @@ public:
         particlePos.setZ( getPeriodicPos( particlePos.getZ(), DirMinusZ, DirPlusZ, (this->MaxBox).getZ(),(this->MinBox).getZ(),DirZ));
     }
 
-    FReal getPeriodicPos(FReal pos, PeriodicCondition periodicPlus, PeriodicCondition periodicMinus,FReal minDir,FReal maxDir,const int dir){
+    FReal getPeriodicPos(FReal pos, PeriodicCondition periodicPlus, PeriodicCondition periodicMinus,FReal maxDir,FReal minDir,const int dir){
         FReal res = pos;
         if( TestPeriodicCondition(dir, periodicPlus) ){
             while(res >= maxDir){

Src/Arranger/FOctreeArranger.hpp

@@ -82,7 +82,7 @@ public:
                 const MortonIndex particuleIndex = tree->getMortonFromPosition(currentPart);
                 if(particuleIndex != currentMortonIndex){
                     //Need to move this one
-                    interface->removeFromLeafAndKeep(particles,idxPart);
+                    interface->removeFromLeafAndKeep(particles,currentPart,idxPart);
                 }
                 else{
                     //Need to increment idx;
@@ -93,7 +93,6 @@ public:
         //Insert back the parts that have been removed
         interface->insertAllParticles(tree);
         //Then, remove the empty leaves
-
         { // Remove empty leaves
             typename OctreeClass::Iterator octreeIterator(tree);
             octreeIterator.gotoBottomLeft();

Src/Kernels/P2P/FP2PLeafInterface.hpp

@@ -40,13 +40,14 @@ public:
 
     void getParticlePosition(ContainerClass* lf,const int idxPart, FPoint& newPos){
         newPos = FPoint(lf->getPositions()[0][idxPart],lf->getPositions()[1][idxPart],lf->getPositions()[2][idxPart]);
+
     }
 
-    void removeFromLeafAndKeep(ContainerClass* lf, const int idxPart){
+    void removeFromLeafAndKeep(ContainerClass* lf, const FPoint& newPos, const int idxPart){
         //Store index
         indexes.push(lf->getIndexes()[idxPart]);
         //Store particles attributes
-        toStoreRemovedParts->push(FPoint(lf->getPositions()[0][idxPart],lf->getPositions()[1][idxPart],lf->getPositions()[2][idxPart]),
+        toStoreRemovedParts->push(newPos,
                                   idxPart,
                                   lf->getPhysicalValues()[idxPart],
                                   lf->getForcesX()[idxPart],lf->getForcesY()[idxPart],lf->getForcesZ()[idxPart]);

Tests/Utils/testOctreeRearrangePeriodic.cpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info".
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+
+#include <iostream>
+
+#include <cstdio>
+#include <cstdlib>
+
+#include "../../Src/Utils/FParameters.hpp"
+#include "../../Src/Utils/FTic.hpp"
+
+#include "../../Src/Containers/FOctree.hpp"
+#include "../../Src/Containers/FVector.hpp"
+
+#include "../../Src/Components/FSimpleLeaf.hpp"
+
+#include "../../Src/Utils/FPoint.hpp"
+
+#include "../../Src/Components/FBasicParticleContainer.hpp"
+#include "../../Src/Components/FBasicCell.hpp"
+#include "../../Src/Kernels/P2P/FP2PLeafInterface.hpp"
+#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
+
+#include "../../Src/Arranger/FOctreeArranger.hpp"
+#include "../../Src/Arranger/FArrangerPeriodic.hpp"
+
+#include "../../Src/Utils/FParameterNames.hpp"
+
+
+
+// Simply create particles and try the kernels
+int main(int argc, char ** argv){
+    FHelpDescribeAndExit(argc, argv,
+                         "Put the particles into a tree, then change the position of some particles and update the tree.\n"
+                         "This method should be used to avoid the tree reconstruction.",
+                         FParameterDefinitions::NbParticles, FParameterDefinitions::OctreeHeight,
+                         FParameterDefinitions::OctreeSubHeight);
+
+    typedef FBasicCell                      CellClass;
+    typedef FP2PParticleContainerIndexed<>      ContainerClass;
+
+    typedef FSimpleLeaf< ContainerClass >                     LeafClass;
+    typedef FOctree< CellClass, ContainerClass , LeafClass >  OctreeClass;
+
+    typedef FArrangerPeriodic<OctreeClass,ContainerClass,FP2PLeafInterface<OctreeClass> > ArrangerClassPeriodic;
+
+    ///////////////////////What we do/////////////////////////////
+    std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
+    //////////////////////////////////////////////////////////////
+
+    const int NbLevels          = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 7);
+    const int SizeSubLevels     = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
+    const int NbPart           = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 2000000);
+
+    FTic counter;
+
+    srand48 ( 1 ); // volontary set seed to constant
+
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+
+    const FReal BoxWidth = 1.0;
+    const FReal BoxCenter = 0.5;
+
+    OctreeClass tree(NbLevels, SizeSubLevels, BoxWidth, FPoint(BoxCenter,BoxCenter,BoxCenter));
+
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+
+    std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
+    std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
+    counter.tic();
+
+
+    {
+
+        FPoint particleToFill;
+        for(int idxPart = 0 ; idxPart < NbPart ; ++idxPart){
+            particleToFill.setPosition(
+                        (BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
+                        (BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)),
+                        (BoxWidth*FReal(drand48())) + (BoxCenter-(BoxWidth/2)));
+            tree.insert(particleToFill,idxPart);
+        }
+    }
+        counter.tac();
+    std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
+
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+    ArrangerClassPeriodic arrange(&tree);
+
+    // In order to test multiple displacement of particles
+    //    for(int ite=0 ; ite<10 ; ++ite){
+    std::cout << "Working on particles ..." << std::endl;
+    counter.tic();
+
+    { // Create new position for each particles
+        OctreeClass::Iterator octreeIterator(&tree);
+        octreeIterator.gotoBottomLeft();
+        do{
+            ContainerClass* particles = octreeIterator.getCurrentListTargets();
+            for(int idxPart = 0; idxPart < particles->getNbParticles() ; ++idxPart){
+                particles->getWPositions()[0][idxPart] = (FReal(drand48()))*BoxWidth*4 + (BoxCenter-(BoxWidth/2));
+                particles->getWPositions()[1][idxPart] = (FReal(drand48()))*BoxWidth*4 + (BoxCenter-(BoxWidth/2));
+                particles->getWPositions()[2][idxPart] = (FReal(drand48()))*BoxWidth*4 + (BoxCenter-(BoxWidth/2));
+            }
+        } while(octreeIterator.moveRight());
+    }
+
+    counter.tac();
+    std::cout << "Done  " << "(@Moving = " << counter.elapsed() << "s)." << std::endl;
+
+
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+
+    std::cout << "Arrange ..." << std::endl;
+    counter.tic();
+
+    //FOctreeArranger<OctreeClass, ContainerClass, TestParticle, Converter<TestParticle> > arrange(&tree);
+    arrange.rearrange();
+
+
+    counter.tac();
+    std::cout << "Done  " << "(@Arrange = " << counter.elapsed() << "s)." << std::endl;
+
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+    //    }
+    std::cout << "Test ..." << std::endl;
+    counter.tic();
+
+    { // Check that each particle has been put into the right leaf
+        long counterPart = 0;
+
+        OctreeClass::Iterator octreeIterator(&tree);
+        octreeIterator.gotoBottomLeft();
+        do{
+            const MortonIndex leafIndex = octreeIterator.getCurrentGlobalIndex();
+
+            ContainerClass* particles = octreeIterator.getCurrentListTargets();
+            for(int idxPart = 0; idxPart < particles->getNbParticles() ; ++idxPart){
+                const FPoint particlePosition( particles->getWPositions()[0][idxPart],
+                                               particles->getWPositions()[1][idxPart],
+                                               particles->getWPositions()[2][idxPart]);
+
+                const MortonIndex particleIndex = tree.getMortonFromPosition( particlePosition );
+                if( leafIndex != particleIndex){
+                    std::cout << "Index problem, should be " << leafIndex <<
+                                 " particleIndex "<< particleIndex << std::endl;
+                }
+
+            }
+
+            counterPart += octreeIterator.getCurrentListTargets()->getNbParticles();
+            if(octreeIterator.getCurrentListTargets()->getNbParticles() == 0){
+                std::cout << "Problem, leaf is empty at index " << leafIndex << std::endl;
+            }
+        } while(octreeIterator.moveRight());
+
+        if( counterPart != NbPart ){
+            std::cout <<"Wrong particles number, should be " << NbPart << " but is " << counterPart << std::endl;
+        }
+    }
+
+    { // Check that each particle has been summed with all other
+        OctreeClass::Iterator octreeIterator(&tree);
+        OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
+
+        const int heightMinusOne = NbLevels - 1;
+        for(int idxLevel = 1 ; idxLevel < heightMinusOne ; ++idxLevel ){
+            // for each cells
+            do{
+                int countChild = 0;
+                CellClass** const child = octreeIterator.getCurrentChild();
+                for(int idxChild = 0 ; idxChild < 8 ; ++idxChild ){
+                    if( child[idxChild] ){
+                        countChild += 1;
+                    }
+                }
+
+                if(countChild == 0){
+                    std::cout << "Problem at level " << idxLevel << " cell has no child " << octreeIterator.getCurrentGlobalIndex() << std::endl;
+                }
+
+            } while(octreeIterator.moveRight());
+
+            avoidGotoLeftIterator.moveDown();
+            octreeIterator = avoidGotoLeftIterator;
+        }
+    }
+
+    counter.tac();
+    std::cout << "Done  " << "(@Test = " << counter.elapsed() << "s)." << std::endl;
+    //////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////////////////
+
+    return 0;
+}