Finished the rotation kernel, BUT need to make an optimized one and remove...

Finished the rotation kernel, BUT need to make an optimized one and remove computational tests. The current version is a copy of the papers formula

Src/Kernels/Rotation/FRotationCell.hpp

+#ifndef FROTATIONCELL_HPP
+#define FROTATIONCELL_HPP
+
+#include "../../Components/FAbstractSerializable.hpp"
+#include "../../Components/FAbstractSendable.hpp"
+#include "../../Utils/FComplexe.hpp"
+#include "../../Utils/FMemUtils.hpp"
+
+#include "../../Extensions/FExtendCellType.hpp"
+
+#include "../../Components/FBasicCell.hpp"
+
+#include "../../Containers/FBufferWriter.hpp"
+#include "../../Containers/FBufferReader.hpp"
+
+/** This class is a cell used for the rotation based kernel
+  * The size of the multipole and local vector are based on a template
+  * User should choose this parameter P carrefuly to match with the
+  * P of the kernel.
+  *
+  * Multipole/Local vectors contain value as:
+  * {0,0}{1,0}{1,1}...{P,P-1}{P,P}
+  * So the size of such vector can be obtained by a suite:
+  * (n+1)*n/2 => (P+2)*(P+1)/2
+  */
+template <int P>
+class FRotationCell : public FBasicCell {
+protected:
+    //< Size of multipole vector
+    static const int MultipoleSize = ((P+2)*(P+1))/2; // Artimethique suite (n+1)*n/2
+    //< Size of local vector
+    static const int LocalSize = ((P+2)*(P+1))/2;     // Artimethique suite (n+1)*n/2
+
+    //< Multipole vector (static memory)
+    FComplexe multipole_exp[MultipoleSize]; //< For multipole extenssion
+    //< Local vector (static memory)
+    FComplexe local_exp[LocalSize];         //< For local extenssion
+
+public:
+    /** Default constructor
+      * Put 0 in vectors
+      */
+    FRotationCell(){
+    }
+
+    /** Copy constructor
+      * Copy the value in the vectors
+      */
+    FRotationCell(const FRotationCell& other){
+        (*this) = other;
+    }
+
+    /** Default destructor */
+    virtual ~FRotationCell(){
+    }
+
+    /** Copy operator
+      * copies only the value in the vectors
+      */
+    FRotationCell& operator=(const FRotationCell& other) {
+        FMemUtils::copyall(multipole_exp, other.multipole_exp, MultipoleSize);
+        FMemUtils::copyall(local_exp, other.local_exp, LocalSize);
+        return *this;
+    }
+
+    /** Get Multipole array */
+    const FComplexe* getMultipole() const {
+        return multipole_exp;
+    }
+    /** Get Local array */
+    const FComplexe* getLocal() const {
+        return local_exp;
+    }
+
+    /** Get Multipole array */
+    FComplexe* getMultipole() {
+        return multipole_exp;
+    }
+    /** Get Local array */
+    FComplexe* getLocal() {
+        return local_exp;
+    }
+
+    ///////////////////////////////////////////////////////
+    // to extend FAbstractSendable
+    ///////////////////////////////////////////////////////
+    void serializeUp(FBufferWriter& buffer) const{
+        buffer.write(multipole_exp, MultipoleSize);
+    }
+    void deserializeUp(FBufferReader& buffer){
+        buffer.fillArray(multipole_exp, MultipoleSize);
+    }
+
+    void serializeDown(FBufferWriter& buffer) const{
+        buffer.write(local_exp, LocalSize);
+    }
+    void deserializeDown(FBufferReader& buffer){
+        buffer.fillArray(local_exp, LocalSize);
+    }
+
+    ///////////////////////////////////////////////////////
+    // to extend Serializable
+    ///////////////////////////////////////////////////////
+    void save(FBufferWriter& buffer) const{
+        FBasicCell::save(buffer);
+        buffer.write(multipole_exp, MultipoleSize);
+        buffer.write(local_exp, LocalSize);
+    }
+    void restore(FBufferReader& buffer){
+        FBasicCell::restore(buffer);
+        buffer.fillArray(multipole_exp, MultipoleSize);
+        buffer.fillArray(local_exp, LocalSize);
+    }
+};
+
+#endif // FROTATIONCELL_HPP

Tests/Kernels/testRotationAlgorithm.cpp

@@ -21,6 +21,7 @@
 #include "../../Src/Core/FFmmAlgorithm.hpp"
 
 #include "../../Src/Kernels/Spherical/FSphericalParticle.hpp"
+#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
 
 #include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
 #include "../../Src/Kernels/Rotation/FRotationCell.hpp"
@@ -53,7 +54,7 @@ public:
 
 
 int main(int argc, char** argv){
-    static const int P = 1;
+    static const int P = 8;
 
     typedef IndexedParticle                ParticleClass;
     typedef FRotationCell<P>               CellClass;
@@ -109,7 +110,7 @@ int main(int argc, char** argv){
     std::cout << "Create kernel ..." << std::endl;
     counter.tic();
 
-    KernelClass kernels(P, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+    KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
 
     counter.tac();
     std::cout << "Done  " << " in " << counter.elapsed() << "s)." << std::endl;
@@ -180,7 +181,7 @@ int main(int argc, char** argv){
                     const ParticleClass& other = particles[leafIter.data().getIndex()];
 
                     potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
-std::cout << leafIter.data().getIndex() << " Direct Potential = " << other.getPotential() << " Fmm Potential = " << leafIter.data().getPotential() << std::endl; // Remove Me
+//std::cout << leafIter.data().getIndex() << " Direct Potential = " << other.getPotential() << " Fmm Potential = " << leafIter.data().getPotential() << std::endl; // Remove Me
                     fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
 
                     fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
@@ -201,6 +202,89 @@ std::cout << leafIter.data().getIndex() << " Direct Potential = " << other.getPo
 
     delete[] particles;
 
+    {
+        typedef FSphericalKernel<ParticleClass, CellClass, ContainerClass>   ShKernelClass;
+        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, ShKernelClass, LeafClass > ShFmmClass;
+
+        FFmaLoader<ParticleClass> loaderSh(filename);
+        OctreeClass treeSh(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+
+        for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+            ParticleClass particle;
+            loaderSh.fillParticle(particle);
+            treeSh.insert(particle);
+        }
+
+        ShKernelClass kernelsSh(P, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+        ShFmmClass algoSh(&treeSh,&kernelsSh);
+        algoSh.execute();
+
+        {
+            const int OctreeHeight = tree.getHeight();
+            typename OctreeClass::Iterator octreeIterator(&treeSh);
+            octreeIterator.gotoBottomLeft();
+
+            typename OctreeClass::Iterator octreeIteratorValide(&tree);
+            octreeIteratorValide.gotoBottomLeft();
+
+            for(int level = OctreeHeight - 1 ; level > 1 ; --level){
+                FMath::FAccurater pole, local;
+
+                do {
+                    for(int idx = 0 ; idx < ((P+2)*(P+1))/2 ; ++idx){
+                        pole.add(octreeIteratorValide.getCurrentCell()->getMultipole()[idx].getReal(),
+                                 octreeIterator.getCurrentCell()->getMultipole()[idx].getReal());
+                        pole.add(octreeIteratorValide.getCurrentCell()->getMultipole()[idx].getImag(),
+                                 octreeIterator.getCurrentCell()->getMultipole()[idx].getImag());
+                    }
+                    for(int idx = 0 ; idx < ((P+2)*(P+1))/2 ; ++idx){
+                        local.add(octreeIteratorValide.getCurrentCell()->getLocal()[idx].getReal(),
+                                  octreeIterator.getCurrentCell()->getLocal()[idx].getReal());
+                        local.add(octreeIteratorValide.getCurrentCell()->getLocal()[idx].getImag(),
+                                  octreeIterator.getCurrentCell()->getLocal()[idx].getImag());
+                    }
+                } while(octreeIterator.moveRight() && octreeIteratorValide.moveRight());
+
+                octreeIterator.moveUp();
+                octreeIterator.gotoLeft();
+
+                octreeIteratorValide.moveUp();
+                octreeIteratorValide.gotoLeft();
+
+                printf("At level %d:", level);
+                printf(">> pole %f, %f\n", pole.getL2Norm(), pole.getInfNorm() );
+                printf(">> local %f, %f\n", local.getL2Norm(), local.getInfNorm() );
+            }
+        }
+        {
+            // Check that each particle has been summed with all other
+            typename OctreeClass::Iterator octreeIterator(&treeSh);
+            octreeIterator.gotoBottomLeft();
+
+            typename OctreeClass::Iterator octreeIteratorValide(&tree);
+            octreeIteratorValide.gotoBottomLeft();
+
+            FMath::FAccurater potential, forces;
+
+            do {
+                typename ContainerClass::BasicIterator iter(*octreeIterator.getCurrentListTargets());
+                typename ContainerClass::BasicIterator iterValide(*octreeIteratorValide.getCurrentListTargets());
+
+                while( iter.hasNotFinished() ){
+
+                   potential.add( iterValide.data().getPotential() , iter.data().getPotential());
+                   forces.add(iterValide.data().getForces().getX(),iter.data().getForces().getX());
+                   forces.add(iterValide.data().getForces().getY(),iter.data().getForces().getY());
+                   forces.add(iterValide.data().getForces().getZ(),iter.data().getForces().getZ());
+
+                    iter.gotoNext();
+                    iterValide.gotoNext();
+                }
+
+            } while(octreeIterator.moveRight() && octreeIteratorValide.moveRight());
+        }
+    }
+
     return 0;
 }
 

UTests/utestRotationDirect.cpp

+
+// ===================================================================================
+// Logiciel initial: ScalFmm Version 0.5
+// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
+// Propriétaires : INRIA.
+// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
+// Initial software: ScalFmm Version 0.5
+// Co-authors: Olivier Coulaud, Bérenger Bramas.
+// Owners: INRIA.
+// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
+// ===================================================================================
+
+#include "../Src/Utils/FGlobal.hpp"
+
+#include "../Src/Containers/FOctree.hpp"
+#include "../Src/Containers/FVector.hpp"
+
+#include "../Src/Kernels/Rotation/FRotationCell.hpp"
+#include "../Src/Kernels/Spherical/FSphericalParticle.hpp"
+
+#include "../Src/Components/FSimpleLeaf.hpp"
+#include "../Src/Kernels/Rotation/FRotationKernel.hpp"
+
+#include "../Src/Files/FFmaBinLoader.hpp"
+#include "../Src/Files/FTreeIO.hpp"
+
+#include "../Src/Core/FFmmAlgorithmThread.hpp"
+#include "../Src/Core/FFmmAlgorithm.hpp"
+
+#include "FUTester.hpp"
+
+/*
+  In this test we compare the spherical fmm results and the direct results.
+  */
+
+/** We need to know the position of the particle in the array */
+class IndexedParticle : public FSphericalParticle {
+    int index;
+public:
+    IndexedParticle(): index(-1){}
+
+    int getIndex() const{
+        return index;
+    }
+    void setIndex( const int inIndex ){
+        index = inIndex;
+    }
+};
+
+/** the test class
+  *
+  */
+class TestRotationDirect : public FUTester<TestRotationDirect> {
+    /** The test method to factorize all the test based on different kernels */
+    template <class ParticleClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
+              class OctreeClass, class FmmClass>
+    void RunTest(const bool isBlasKernel){
+        // Warning in make test the exec dir it Build/UTests
+        // Load particles
+        FFmaBinLoader<ParticleClass> loader("../../Data/utestSphericalDirect.bin.fma");
+        if(!loader.isOpen()){
+            Print("Cannot open particles file.");
+            uassert(false);
+            return;
+        }
+        Print("Number of particles:");
+        Print(loader.getNumberOfParticles());
+
+        const int NbLevels      = 4;
+        const int SizeSubLevels = 2;
+
+        // Create octree
+        OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+        ParticleClass* const particles = new ParticleClass[loader.getNumberOfParticles()];
+        for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+            loader.fillParticle(particles[idxPart]);
+            particles[idxPart].setIndex( idxPart );
+            tree.insert(particles[idxPart]);
+        }
+
+
+        // Run FMM
+        Print("Fmm...");
+        //KernelClass kernels(NbLevels,loader.getBoxWidth());
+        KernelClass kernels(NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
+        FmmClass algo(&tree,&kernels);
+        algo.execute();
+
+        // Run direct computation
+        Print("Direct...");
+        for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
+            for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
+                kernels.particlesMutualInteraction(&particles[idxTarget], &particles[idxOther]);
+            }
+        }
+
+        // Compare
+        Print("Compute Diff...");
+        FMath::FAccurater potentialDiff;
+        FMath::FAccurater fx, fy, fz;
+        { // Check that each particle has been summed with all other
+            typename OctreeClass::Iterator octreeIterator(&tree);
+            octreeIterator.gotoBottomLeft();
+
+            do{
+                typename ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());
+
+                while( leafIter.hasNotFinished() ){
+                    const ParticleClass& other = particles[leafIter.data().getIndex()];
+
+                    potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
+
+                    fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
+
+                    fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
+
+                    fz.add(other.getForces().getZ(),leafIter.data().getForces().getZ());
+
+                    leafIter.gotoNext();
+                }
+            } while(octreeIterator.moveRight());
+        }
+
+        delete[] particles;
+
+        // Print for information
+        Print("Potential diff is = ");
+        Print(potentialDiff.getL2Norm());
+        Print(potentialDiff.getInfNorm());
+        Print("Fx diff is = ");
+        Print(fx.getL2Norm());
+        Print(fx.getInfNorm());
+        Print("Fy diff is = ");
+        Print(fy.getL2Norm());
+        Print(fy.getInfNorm());
+        Print("Fz diff is = ");
+        Print(fz.getL2Norm());
+        Print(fz.getInfNorm());
+
+        // Assert
+        const FReal MaximumDiff = FReal(0.0001);
+        uassert(potentialDiff.getL2Norm() < MaximumDiff);
+        uassert(potentialDiff.getInfNorm() < MaximumDiff);
+        uassert(fx.getL2Norm()  < MaximumDiff);
+        uassert(fx.getInfNorm() < MaximumDiff);
+        uassert(fy.getL2Norm()  < MaximumDiff);
+        uassert(fy.getInfNorm() < MaximumDiff);
+        uassert(fz.getL2Norm()  < MaximumDiff);
+        uassert(fz.getInfNorm() < MaximumDiff);
+    }
+
+    /** If memstas is running print the memory used */
+    void PostTest() {
+        if( FMemStats::controler.isUsed() ){
+            std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
+            std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
+            std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
+        }
+    }
+
+    ///////////////////////////////////////////////////////////
+    // The tests!
+    ///////////////////////////////////////////////////////////
+
+    static const int P = 9;
+
+    /** Rotation */
+    void TestRotation(){
+        typedef IndexedParticle         ParticleClass;
+        typedef FRotationCell<P>            CellClass;
+        typedef FVector<ParticleClass>  ContainerClass;
+
+        typedef FRotationKernel<ParticleClass, CellClass, ContainerClass, P >          KernelClass;
+
+        typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
+        typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;
+
+        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
+
+        RunTest<ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass,
+                OctreeClass, FmmClass>(false);
+    }
+
+    ///////////////////////////////////////////////////////////
+    // Set the tests!
+    ///////////////////////////////////////////////////////////
+
+    /** set test */
+    void SetTests(){
+        AddTest(&TestRotationDirect::TestRotation,"Test Rotation Kernel");
+    }
+};
+
+
+// You must do this
+TestClass(TestRotationDirect)
+
+
+