testFmmDemonstration.cpp 9.25 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 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.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12
//
13 14 15
// 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
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19
// ===================================================================================
20 21 22

#include <iostream>

23 24
#include <cstdio>
#include <cstdlib>
25

26 27 28
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FGlobal.hpp"
29

30 31
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
32

33
#include "../../Src/Components/FSimpleLeaf.hpp"
34

35
#include "../../Src/Utils/FPoint.hpp"
36

37 38
#include "../../Src/Components/FTestCell.hpp"
#include "../../Src/Components/FTestKernels.hpp"
39

40 41
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
42 43


44
#include "../../Src/Components/FAbstractParticleContainer.hpp"
45
#include "../../Src/Components/FBasicKernels.hpp"
46

47
#include "../../Src/Files/FRandomLoader.hpp"
48

49 50
#include "../../Src/Utils/FParameterNames.hpp"

51 52 53 54
// My cell is actually a basic cell => minimum of data
class MyCell : public FBasicCell {
};

55 56
// My fack Container is simply saving the indexes of each
// particles (nothing more!)
57 58
template <class FReal>
class MyContainer : public FAbstractParticleContainer<FReal> {
59
    FVector<FSize> indexes;
60
public:
61
    template<typename... Args>
62
    void push(const FPoint<FReal>& /*inParticlePosition*/, const FSize newParticleIndex, Args ... /*args*/){
63
        indexes.push(newParticleIndex);
64
    }
65

66
    FSize getSize() const {
67
        return indexes.getSize();
68
    }
69

70
    const FVector<FSize>& getIndexes() const{
71
        return indexes;
72 73 74
    }
};

75 76
// My leaf process the particles and save only
// those where keepIt is true (during the push method)
77 78 79
template <class FReal>
class MyLeaf : public FAbstractLeaf< FReal, MyContainer<FReal> > {
    MyContainer<FReal> particles;
80

81
public:
82
    template<typename... Args>
83
    void push(const FPoint<FReal>& inParticlePosition, const bool keepIt, Args ... args){
84
        if(keepIt) particles.push(inParticlePosition, args...);
85
    }
86
    MyContainer<FReal>* getSrc(){
87
        return &particles;
88
    }
89
    MyContainer<FReal>* getTargets(){
90
        return &particles;
91 92 93 94 95
    }
};

// My kernel actually does nothing except showing how to retreive data from an
// array from the indexes vector giving by the leaf in the P2M
96 97 98
template< class CellClass, class ContainerClass>
class MyKernel : public FAbstractKernels<CellClass,ContainerClass>{
    MortonIndex* indexForEachParticle;
99
public:
100
    MyKernel(const FSize inNbParticles): indexForEachParticle(new MortonIndex[inNbParticles]) {
101
        memset(indexForEachParticle,0,sizeof(MortonIndex)*inNbParticles);
102 103
    }

104 105 106 107
    ~MyKernel(){
        delete[] indexForEachParticle;
    }

108
    void P2M(CellClass* const cell, const ContainerClass* const particles)  override  {
109
        for(FSize idxPart = 0 ; idxPart < particles->getSize() ; ++idxPart){
110 111
            // save the current morton index for each particles
            indexForEachParticle[ particles->getIndexes()[idxPart] ] = cell->getMortonIndex();
112 113 114
        }
    }

115
    void M2M(CellClass* const FRestrict , const CellClass*const FRestrict *const FRestrict , const int )  override  {
116 117
    }

118
    void M2L(CellClass* const FRestrict , const CellClass* [], const int [], const int , const int )  override  {
119 120
    }

121
    void L2L(const CellClass* const FRestrict , CellClass* FRestrict *const FRestrict  , const int ) override  {
122 123
    }

124
    void L2P(const CellClass* const , ContainerClass* const ) override {
125 126
    }

berenger-bramas's avatar
berenger-bramas committed
127
    void P2P(const FTreeCoordinate& ,
128 129
                     ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
                     ContainerClass* const [], const int [], const int ) override {
berenger-bramas's avatar
berenger-bramas committed
130

131 132
    }

133 134 135 136 137

    void P2POuter(const FTreeCoordinate& ,
                     ContainerClass* const FRestrict ,
                     ContainerClass* const [],
                    const int [], const int ) override {
berenger-bramas's avatar
berenger-bramas committed
138

139
    }
140

141 142 143 144
};


int main(int argc, char ** argv){
145 146 147 148 149
    FHelpDescribeAndExit(argc, argv,
                         "Explains how to use ScalFMM (only the code is interesting).",
                         FParameterDefinitions::OctreeHeight, FParameterDefinitions::OctreeSubHeight,
                         FParameterDefinitions::NbParticles);

150
    typedef double FReal;
151
    // Custom data structure here
152
    typedef MyCell            CellClass;
153 154
    typedef MyContainer<FReal>       ContainerClass;
    typedef MyLeaf<FReal>            LeafClass;
155

156
    // Standard things here
157
    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
158 159
    typedef MyKernel< CellClass, ContainerClass >         KernelClass;
    typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass >     FmmClass;
160 161 162 163 164 165 166 167 168

    //////////////////////////////////////////////////////////////
    ///////////////////////What we do/////////////////////////////

    std::cout << ">> This executable has to be used to demonstrate the use of scalfmm.\n";

    //////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////

169 170
    const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
    const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
171
    const FSize NbPart = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, FSize(20));
172 173 174 175 176
    FTic counter;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

177
    FRandomLoader<FReal> loader(NbPart, 1, FPoint<FReal>(0.5,0.5,0.5), 1);
178
    OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
179 180 181 182 183 184 185 186

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
    std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
    counter.tic();

187
    FPoint<FReal>*const realsParticlesPositions = new FPoint<FReal>[NbPart];
188
    {
189
        FPoint<FReal> particlePosition;
190
        bool keepIt;
191
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
192
            // get a random position
193 194 195
            loader.fillParticle(&particlePosition);
            // let say we remove 1/5 particles
            keepIt = (idxPart%5);
196
            // insert in the tree
197 198 199
            tree.insert(particlePosition, keepIt, idxPart);
            // save the position
            realsParticlesPositions[idxPart] = particlePosition;
200 201 202 203 204 205 206 207 208 209 210 211
        }
    }

    counter.tac();
    std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Which particles are in wich leafs ..." << std::endl;
    counter.tic();

212
    OctreeClass::Iterator octreeIterator(&tree);
213 214
    octreeIterator.gotoBottomLeft();
    do{
215
        FVector<FSize>::ConstBasicIterator iter(octreeIterator.getCurrentListTargets()->getIndexes());
216 217 218
        const MortonIndex indexAtThisLeaf = octreeIterator.getCurrentGlobalIndex();

        while( iter.hasNotFinished() ){
219 220 221
            std::cout << "Particles with index " << iter.data() <<
                         " has a morton index of " << indexAtThisLeaf << std::endl;
            std::cout << " it real position was " << realsParticlesPositions[iter.data()] << std::endl;
222 223 224 225 226 227 228 229 230 231 232 233 234
            iter.gotoNext();
        }
    } while(octreeIterator.moveRight());

    counter.tac();
    std::cout << "Done  " << "(@Counting = " << counter.elapsed() << "s)." << std::endl;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

    std::cout << "Working on particles ..." << std::endl;
    counter.tic();

235
    KernelClass kernels(NbPart);
236 237 238 239 240 241 242 243 244
    FmmClass algo(&tree,&kernels);
    algo.execute();

    counter.tac();
    std::cout << "Done  " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;

    //////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////

245
    delete [] realsParticlesPositions;
246 247 248 249 250

    return 0;
}


251