FUserKernelEngine.hpp 34.6 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5 6 7 8 9 10

/**
 * @file This file contains a class that inherits from FScalFMMEngine,
 * and will implement the API functions for a user defined kernel.
 */
#ifndef FUSERKERNELENGINE_HPP
#define FUSERKERNELENGINE_HPP

#include "FScalFMMEngine.hpp"
11
#include "FUserLeafContainer.hpp"
12
#include <vector>
13 14 15 16

/**
 * @brief CoreCell : Cell used to store User datas
 */
17
class CoreCell : public FBasicCell, public FExtendCellType {
18 19 20
    // Mutable in order to work with the API
    mutable void* userData;

21
protected:
22 23 24
    //Static members to be initialised before octree creation
    static Scalfmm_Cell_Descriptor user_cell_descriptor;

25
public:
26 27 28 29 30 31 32 33 34 35 36 37
    static void Init(Scalfmm_Cell_Descriptor cell_descriptor){
        user_cell_descriptor=cell_descriptor;
    }

    static Callback_init_cell GetInit(){
        return user_cell_descriptor.user_init_cell;
    }

    static Callback_free_cell GetFree(){
        return user_cell_descriptor.user_free_cell;
    }

38 39 40
    CoreCell() : userData(nullptr) {
    }

41
    //We free the cells here
42
    ~CoreCell(){
43 44 45
        if(userData){
            this->user_cell_descriptor.user_free_cell(userData);
        }
46
    }
47

48 49 50 51 52 53 54 55 56 57 58 59 60
    /**
     * @brief setContainer store the ptr to the user data inside our
     * struct
     */
    void setContainer(void* inContainer) const {
        userData = inContainer;
    }

    /**
     * @brief getContainer : return the user datas (in order to give
     * it back to the user defined kernel function)
     */
    void* getContainer() const {
61 62 63 64 65 66 67
        if(userData){
            return userData;
        }
        else{
            //std::cout << "UserDatas not initialised\n";
            return nullptr; //pareil que userdata, finalement
        }
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
    }
};


/**
 * This class simply call the function pointers from Scalfmm_Kernel_Descriptor.
 * If not pointer is set the calls are skipped.
 * The userData is given at any calls.
 */
template< class CellClass, class ContainerClass>
class CoreKernel : public FAbstractKernels<CellClass,ContainerClass> {
    Scalfmm_Kernel_Descriptor kernel;
    void* userData;

public:
    CoreKernel(Scalfmm_Kernel_Descriptor inKernel, void* inUserData) : kernel(inKernel) , userData(inUserData){
    }

    /** Default destructor */
    virtual ~CoreKernel(){
88

89 90 91 92
    }

    /** Do nothing */
    virtual void P2M(CellClass* const cell, const ContainerClass* const container) {
93
        if(kernel.p2m) kernel.p2m(cell->getContainer(), container->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData);
94 95 96 97
    }

    /** Do nothing */
    virtual void M2M(CellClass* const FRestrict cell, const CellClass*const FRestrict *const FRestrict children, const int level) {
98 99
        if(kernel.m2m_full){
            std::vector<void *> userCellArray;
100 101 102 103 104 105
            for(int idx=0 ;idx<8 ; ++idx){
                if( children[idx] ){
                    userCellArray.push_back(children[idx]->getContainer());
                }else{
                    userCellArray.push_back(nullptr);
                }
106
            }
107
            kernel.m2m_full(level, cell->getContainer(), userCellArray.data(), userData);
108 109 110 111 112 113
        }else{
            if(kernel.m2m){
                for(int idx = 0 ; idx < 8 ; ++idx){
                    if( children[idx] ){
                        kernel.m2m(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
                    }
114 115 116 117 118 119
                }
            }
        }
    }

    /** Do nothing */
120 121
    virtual void M2L(CellClass* const FRestrict cell, const CellClass* distanNeighbors[],
                     const int neighborPositions[], const int size,const int level) {
122 123
        if(kernel.m2l_full){//all 343 interactions will be computed directly
            //First, copy the fmm cell inside an array of user cells
124 125 126 127
            std::vector<void *> userCellArray;
            userCellArray.resize(size);
            for(int i=0 ; i<size ; ++i){
                userCellArray[i] = distanNeighbors[i]->getContainer();
128
            }
129 130
            kernel.m2l_full(level,cell->getContainer(),neighborPositions,size,userCellArray.data(),userData);
            FAssertLF("m2l_full temporary disabled ...\n");
131
        }else{
132
            if(kernel.m2l){
133 134 135
                for(int idx = 0 ; idx < size ; ++idx){
                    const int idxNeigh = neighborPositions[idx];
                    kernel.m2l(level, cell->getContainer(), idxNeigh, distanNeighbors[idx]->getContainer(), userData);
136 137 138 139 140 141
                }
            }
        }
    }

    /** Do nothing */
142 143
    virtual void L2L(const CellClass* const FRestrict cell,
                     CellClass* FRestrict *const FRestrict children, const int level) {
144 145
        if(kernel.l2l_full){
            std::vector<void *> userCellArray;
146
            std::cout<<"Children array is good ?\t"<<children<<"\n";
147
            for(int i=0 ;i<8 ; ++i){
148
                if(children[i]){
149
                    std::cout<<"Children["<<i<<"] is good ?\t"<<children[i]<<"\n";
150 151 152 153
                    userCellArray.push_back(children[i]->getContainer());
                }else{
                    userCellArray.push_back(nullptr);
                }
154
            }
155
            kernel.l2l_full(level, cell->getContainer(), userCellArray.data(), userData);
156 157 158 159 160 161
        }else{
            if(kernel.l2l){
                for(int idx = 0 ; idx < 8 ; ++idx){
                    if( children[idx] ){
                        kernel.l2l(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
                    }
162 163 164 165 166 167
                }
            }
        }
    }

    virtual void L2P(const CellClass* const cell, ContainerClass* const container){
168
        //        std::cout << "L2P with "<< container->getNbParticles() <<" - over "<< cell<< " and "<<cell->getContainer() <<" Indexes : ["<<container->getIndexes()[0] <<"]\n";
169
        if(kernel.l2p) kernel.l2p(cell->getContainer(),container->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData);
170 171
    }

172 173 174 175
    virtual void P2POuter(const FTreeCoordinate& inLeafPosition,
                          ContainerClass* const FRestrict targets,
                          ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
                          const int size){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
176 177 178 179
        // for(int idx = 0 ; idx < size ; ++idx){
        //     kernel.p2p(targets->getNbParticles(),targets.getIndexes().data(),
        //                directNeighborsParticles[idx]->getNbParticles(),
        // }
180
    }
181

182 183 184
    virtual void P2P(const FTreeCoordinate& inLeafPosition,
                     ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                     ContainerClass* const neighbors[],const int sourcePosition[], const int size){
185 186 187 188
        if(kernel.p2pinner){
            kernel.p2pinner(targets->getContainer(),targets->getNbParticles(), targets->getIndexes().data(), userData);
            //std::cout<<"Inner used\n";
        }
189 190
        if(kernel.p2p_full){
            //Create the arrays of size and indexes
191 192 193 194 195 196 197 198
            FSize * nbPartPerNeighbors = new FSize[size];
            const FSize ** indicesPerNeighbors = new const FSize*[size];
            //create artificial array of void * :
            void ** arrayOfUserContainer = new void *[size];
            for(int idx=0 ; idx<size ; ++idx){
                nbPartPerNeighbors[idx] = neighbors[idx]->getNbParticles();
                indicesPerNeighbors[idx] = neighbors[idx]->getIndexes().data();
                arrayOfUserContainer[idx] = neighbors[idx]->getContainer();
199
            }
200
            kernel.p2p_full(targets->getContainer(),targets->getNbParticles(),targets->getIndexes().data(),arrayOfUserContainer,indicesPerNeighbors,nbPartPerNeighbors,sourcePosition,size,userData);
201 202 203
            delete [] nbPartPerNeighbors;
            delete [] indicesPerNeighbors;
            delete [] arrayOfUserContainer;
204
            //            std::cout<<"P2PFull used\n";
205
        }
206
        if(kernel.p2p_sym){
207
            for(int idx = 0 ; ((idx < size) && (sourcePosition[idx] < 14)) ; ++idx){
208 209
                kernel.p2p_sym(targets->getContainer(),targets->getNbParticles(), targets->getIndexes().data(),
                               neighbors[idx]->getContainer(),neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
210
            }
211
            //            std::cout<<"P2PSym used\n";
212 213 214
        }
        else{
            if(kernel.p2p){
215
                for(int idx = 0 ; idx < size ; ++idx){
216 217
                    kernel.p2p(targets->getContainer(),targets->getNbParticles(), targets->getIndexes().data(),
                               neighbors[idx]->getContainer(),neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
218
                }
219
                //                std::cout<<"P2P used\n";
220 221 222 223
            }
        }
    }

224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
    virtual void P2PRemote(const FTreeCoordinate& inLeafPosition,
                           ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                           const ContainerClass* const neighbors[],const int sourcePosition[], const int size){
        if(kernel.p2p_remote){
            //Create the arrays of size and indexes
            FSize * nbPartPerNeighbors = new FSize[size];
            const FSize ** indicesPerNeighbors = new const FSize*[size];
            //create artificial array of void * :
            void ** arrayOfUserContainer = new void *[size];
            for(int idx=0 ; idx<size ; ++idx){
                nbPartPerNeighbors[idx] = neighbors[idx]->getNbParticles();
                indicesPerNeighbors[idx] = neighbors[idx]->getIndexes().data();
                arrayOfUserContainer[idx] = neighbors[idx]->getContainer();
            }
            kernel.p2p_remote(targets->getContainer(),targets->getNbParticles(),targets->getIndexes().data(),arrayOfUserContainer,
                              indicesPerNeighbors,nbPartPerNeighbors,sourcePosition,size,userData);
            delete [] nbPartPerNeighbors;
            delete [] indicesPerNeighbors;
            delete [] arrayOfUserContainer;
        }else{
            std::cout<<"Warning, MPI used but no P2P Remote set\n";
        }
246 247
    }

248 249 250 251
    //Getter
    void * getUserKernelDatas(){
        return userData;
    }
252 253 254 255
    //Getter
    Scalfmm_Kernel_Descriptor getKernelFct() const {
        return kernel;
    }
256

257 258 259 260 261 262 263
    void M2L_Extended(CellClass * src, CellClass * tgt, const FTreeCoordinate transfer, const int level){
        if(kernel.m2l_ext){
            int array[3] = {transfer.getX(),transfer.getY(),transfer.getZ()};
            kernel.m2l_ext(level,tgt->getContainer(),src->getContainer(),array,userData);
        }
    }

264 265
};

266 267
template<class FReal,class LeafClass>
class FUserKernelEngine : public FScalFMMEngine<FReal>{
268 269

private:
270
    //Typedefs
271
    using ContainerClass = FUserLeafContainer<FReal>;
272 273

    //Typedefs :
274 275
    using OctreeClass = FOctree<FReal,CoreCell,ContainerClass,LeafClass>;
    using CoreKernelClass =  CoreKernel<CoreCell,ContainerClass>;
276

277
    //For arranger classes
278 279 280 281

    //Attributes
    OctreeClass * octree;
    CoreKernelClass * kernel;
282
    int upperLimit;
283 284 285
    // ArrangerClass * arranger;
    // ArrangerClassTyped * arrangerTyped;
    // ArrangerClassPeriodic * arrangerPeriodic;
286

287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
protected:
    int treeHeight;
    FPoint<FReal> boxCenter;
    FPoint<FReal> boxCorner;
    FReal boxWidth;
    FReal boxWidthAtLeafLevel;


    FPoint<FReal> getBoxCenter() const {
        return boxCenter;
    }
    FPoint<FReal> getBoxCorner() const {
        return boxCorner;
    }
    FReal getBoxWidth() const {
        return boxWidth;
    }
    FReal getBoxWidthAtLeafLevel() const {
        return boxWidthAtLeafLevel;
    }
    int getTreeHeight() const {
        return treeHeight;
    }
    OctreeClass* getTree() const {
        return octree;
    }

314

315
public:
316
    FUserKernelEngine(/*int TreeHeight, double BoxWidth , double * BoxCenter, */scalfmm_kernel_type KernelType, scalfmm_algorithm algo) :
317
        octree(nullptr), kernel(nullptr), upperLimit(2),treeHeight(0), boxCenter(0,0,0), boxCorner(0,0,0), boxWidth(0) /*,arranger(nullptr)*/ {
318
        FScalFMMEngine<FReal>::kernelType = KernelType;
319 320 321 322
        FScalFMMEngine<FReal>::Algorithm = algo;
    }

    FUserKernelEngine() : octree(nullptr), kernel(nullptr), upperLimit(2), treeHeight(0) {
323 324
    }

325

326 327
    ~FUserKernelEngine(){
        delete octree;
328
        octree=nullptr;
329 330 331 332
        // if(arranger){
        //     delete arranger;
        //     arranger=nullptr;
        // }
333 334
        if(kernel){
            delete kernel;
335
            kernel=nullptr;
336 337
        }
    }
338

339 340 341 342
    CoreKernelClass * getKernelPtr() {
        return kernel;
    }

343
    virtual void user_kernel_config( Scalfmm_Kernel_Descriptor userKernel, void * userDatas){
344 345 346
        if(!kernel){
            kernel = new CoreKernelClass(userKernel,userDatas);
        }
347 348
    }

349 350 351
    virtual void build_tree(int TreeHeight,double BoxWidth,double* BoxCenter,
                            Scalfmm_Cell_Descriptor user_cell_descriptor,
                            Scalfmm_Leaf_Descriptor user_leaf_descriptor){
352
        CoreCell::Init(user_cell_descriptor);
353
        ContainerClass::Init(user_leaf_descriptor);
354
        this->treeHeight = TreeHeight;
355 356 357 358 359 360
        this->boxCenter = FPoint<FReal>(BoxCenter[0],BoxCenter[1],BoxCenter[2]);
        this->boxWidth = BoxWidth;
        boxCorner.setX(boxCenter.getX() - boxWidth/2.0);
        boxCorner.setY(boxCenter.getY() - boxWidth/2.0);
        boxCorner.setZ(boxCenter.getZ() - boxWidth/2.0);
        this->boxWidthAtLeafLevel = BoxWidth/(2<<TreeHeight);
361
        printf("Tree Height : %d \n",TreeHeight);
362
        this->octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint<FReal>(BoxCenter));
363
    }
364

365
    virtual void apply_on_cell(Callback_apply_on_cell function){
366 367 368 369 370 371
        double boxwidth = octree->getBoxWidth();
        //apply user function reset on each user's cell
        octree->forEachCellWithLevel([&](CoreCell * currCell,const int currLevel){
                if(currCell->getContainer()){
                    FTreeCoordinate currCoord = currCell->getCoordinate();
                    int arrayCoord[3] = {currCoord.getX(),currCoord.getY(),currCoord.getZ()};
372
                    MortonIndex    currMorton = currCoord.getMortonIndex();
373
                    double position[3];
374 375 376
                    position[0] = boxCorner.getX() + currCoord.getX()*boxwidth/double(1<<currLevel);
                    position[1] = boxCorner.getY() + currCoord.getY()*boxwidth/double(1<<currLevel);
                    position[2] = boxCorner.getZ() + currCoord.getZ()*boxwidth/double(1<<currLevel);
377
                    function(currLevel,currMorton,arrayCoord,position,currCell->getContainer(),kernel->getUserKernelDatas());
378 379 380 381 382 383 384 385 386 387 388 389 390 391
                }
            });
    }


    void tree_insert_particles( int NbPositions, double * X, double * Y, double * Z, PartType type){
        if(type == BOTH){
            for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
                octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),idPart);
            }
            FScalFMMEngine<FReal>::nbPart += NbPositions;
        }else{
            if(type==SOURCE){
                for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
392
                    octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeSource,idPart);
393 394 395 396
                }
                FScalFMMEngine<FReal>::nbPart += NbPositions;
            }else{
                for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
397
                    octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),FParticleType::FParticleTypeTarget,idPart);
398 399 400
                }
                FScalFMMEngine<FReal>::nbPart += NbPositions;
            }
401
        }
402

403
        this->init_cell();
404 405
    }

406 407 408 409 410 411 412 413 414
    void tree_insert_particles_xyz( int NbPositions, double * XYZ, PartType type){
        if(type == BOTH){
            for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
                octree->insert(FPoint<FReal>(&XYZ[3*idPart]),idPart);
            }
            FScalFMMEngine<FReal>::nbPart += NbPositions;
        }else{
            if(type==SOURCE){
                for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
415
                    octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeSource,idPart);
416 417 418 419
                }
                FScalFMMEngine<FReal>::nbPart += NbPositions;
            }else{
                for(FSize idPart = 0; idPart<NbPositions ; ++idPart){
420
                    octree->insert(FPoint<FReal>(&XYZ[3*idPart]),FParticleType::FParticleTypeTarget,idPart);
421 422 423
                }
                FScalFMMEngine<FReal>::nbPart += NbPositions;
            }
424
        }
425
        this->init_cell();
426 427
    }

428 429 430
    /**
     * To retrieve the positions, in order to move the parts
     */
431 432
    void get_positions_xyz(int NbPositions, double * positionsToFill, PartType type){
        FScalFMMEngine<FReal>::template generic_get_positions_xyz<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,positionsToFill,type);
433 434
    }

435 436
    void get_positions_xyz_npart(int NbPositions, int * idxOfParticles, double * positionsToFill,PartType type){
        FScalFMMEngine<FReal>::template generic_get_positions_xyz_npart<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,idxOfParticles,positionsToFill,type);
437 438
    }

439
    void get_positions(int NbPositions, double *X, double *Y , double *Z, PartType type){
440
        FScalFMMEngine<FReal>::template generic_get_positions<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,X,Y,Z,type);
441 442
    }

443 444
    void get_positions_npart(int NbPositions, int * idxOfParticles,double * X, double * Y , double * Z,PartType type){
        FScalFMMEngine<FReal>::template generic_get_positions_npart<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,idxOfParticles,X,Y,Z,type);
445 446 447 448 449 450
    }



    //Arranger parts : following function provide a way to move parts
    //inside the tree
451 452
    void add_to_positions_xyz(int NbPositions,double * updatedXYZ,PartType type){
        FScalFMMEngine<FReal>::template generic_add_to_positions_xyz<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,updatedXYZ,type);
453 454
    }

455 456
    void add_to_positions(int NbPositions,double * X, double * Y , double * Z,PartType type){
        FScalFMMEngine<FReal>::template generic_add_to_positions<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,X,Y,Z,type);
457 458
    }

459 460 461
    void set_positions_xyz(int NbPositions, FReal * updatedXYZ, PartType type){
        FScalFMMEngine<FReal>::template generic_set_positions_xyz<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,updatedXYZ,type);
    }
462

463 464
    void set_positions(int NbPositions, FReal * X,FReal * Y,FReal * Z, PartType type){
        FScalFMMEngine<FReal>::template generic_set_positions<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,X,Y,Z,type);
465 466
    }

467 468 469 470 471
    void set_positions_xyz_npart(int NbPositions, int* idxOfParticles, FReal * updatedXYZ, PartType type){
        FScalFMMEngine<FReal>::template generic_set_positions_xyz_npart<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,idxOfParticles,updatedXYZ,type);
    }
    void set_positions_npart(int NbPositions, int* idxOfParticles, FReal * X, FReal * Y , FReal * Z, PartType type){
        FScalFMMEngine<FReal>::template generic_set_positions_npart<ContainerClass,LeafClass,CoreCell>(octree,NbPositions,idxOfParticles,X,Y,Z,type);
472 473
    }

474
    virtual void apply_on_each_leaf(Callback_apply_on_leaf function){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
475 476 477
        generic_apply_on_each_leaf<ContainerClass,CoreCell>(octree,kernel->getUserKernelDatas(),function);
    }

478 479 480
    template<class ContainerClass,class CellClass>
    void generic_apply_on_each_leaf(FOctree<FReal,CellClass,ContainerClass,LeafClass>* octreeIn,
                                    void * kernelUserData,
481
                                    Callback_apply_on_leaf function){
482 483 484 485 486 487 488 489 490 491 492 493
        if(octreeIn){
            octreeIn->forEachCellLeaf([&](CoreCell * currCell, LeafClass * leaf){
                    int lvl = octreeIn->getHeight();
                    MortonIndex currMorton = currCell->getMortonIndex();
                    //Computation of the Center from Coordinate
                    FTreeCoordinate treeCoord = currCell->getCoordinate();
                    double boxWidthLeafLevel = octreeIn->getBoxWidth() / (2 << lvl);
                    FPoint<double> absolutCoord = FPoint<double>(treeCoord.getX()*boxWidthLeafLevel,
                                                                 treeCoord.getY()*boxWidthLeafLevel,
                                                                 treeCoord.getZ()*boxWidthLeafLevel);
                    FPoint<double> leafCenter = absolutCoord + (octreeIn->getBoxCenter()-octreeIn->getBoxWidth()) + boxWidthLeafLevel/2;
                    function(lvl,leaf->getSrc()->getNbParticles(),leaf->getSrc()->getIndexes().data(),currMorton,leafCenter.getDataValue(),
494
                             currCell->getContainer(),leaf->getSrc()->getContainer(),kernelUserData);
495 496 497 498 499
                });
        }else{
            std::cout << "Need to Build the tree and insert the parts First\n" << std::endl;
        }
    }
500

501 502 503
    /*
     * Call the user allocator on userDatas member field of each cell
     */
504
    virtual void init_cell(){
505 506 507 508 509 510 511
        void * generic_ptr = nullptr;
        if(kernel){
            generic_ptr = kernel->getUserKernelDatas();
        }
        else{
            std::cout <<"Warning, no user kernel data set, need to call kernel config first"<< std::endl;
        }
512 513 514 515
        double boxwidth = octree->getBoxWidth();
        //apply user function on each cell
        octree->forEachCellWithLevel([&](CoreCell * currCell,const int currLevel){
                if(!(currCell->getContainer())){
516 517
                    FTreeCoordinate currCoord = currCell->getCoordinate();
                    int arrayCoord[3] = {currCoord.getX(),currCoord.getY(),currCoord.getZ()};
518
                    MortonIndex    currMorton = currCoord.getMortonIndex();
519
                    double position[3];
520 521 522
                    position[0] = boxCorner.getX() + currCoord.getX()*boxwidth/double(1<<currLevel);
                    position[1] = boxCorner.getY() + currCoord.getY()*boxwidth/double(1<<currLevel);
                    position[2] = boxCorner.getZ() + currCoord.getZ()*boxwidth/double(1<<currLevel);
523
                    currCell->setContainer(CoreCell::GetInit()(currLevel,currMorton,arrayCoord,position,generic_ptr));
524 525
                }
            });
526
        //Then init leaves
527 528 529 530 531 532 533 534 535 536 537 538 539 540 541
        if(ContainerClass::GetInitLeaf()){
            octree->forEachCellLeaf([&](CoreCell * currCell, LeafClass * leaf){
                    FTreeCoordinate currCoord = currCell->getCoordinate();
                    int currLevel = octree->getHeight()-1;
                    MortonIndex    currMorton = currCoord.getMortonIndex();
                    double position[3];
                    position[0] = boxCorner.getX() + currCoord.getX()*boxwidth/double(1<<currLevel);
                    position[1] = boxCorner.getY() + currCoord.getY()*boxwidth/double(1<<currLevel);
                    position[2] = boxCorner.getZ() + currCoord.getZ()*boxwidth/double(1<<currLevel);
                    leaf->getSrc()->setContainer(ContainerClass::GetInitLeaf()(currLevel,
                                                                               leaf->getSrc()->getNbParticles(),
                                                                               leaf->getSrc()->getIndexes().data(), currMorton,
                                                                               position, currCell->getContainer(), this->kernel->getUserKernelDatas()));
                });
        }
542 543
    }

544

545
    void free_cell(Callback_free_cell user_cell_deallocator, Callback_free_leaf free_leaf){
546 547 548 549 550 551
        if(free_leaf){
            octree->forEachCellLeaf([&](CoreCell * currCell, LeafClass * leaf){
                    free_leaf(currCell->getContainer(),leaf->getSrc()->getNbParticles(), leaf->getSrc()->getIndexes().data(),
                              leaf->getSrc()->getContainer(),this->kernel->getUserKernelDatas());
                });
        }
552
        octree->forEachCell([&](CoreCell * currCell){
553 554
                if(currCell->getContainer()){
                    user_cell_deallocator(currCell->getContainer());
555
                    currCell->setContainer(nullptr);
556
                }
557 558 559
            });
    }

560 561 562 563 564 565 566 567 568 569
    void set_upper_limit(int inUpperLimit){
        upperLimit = inUpperLimit;
    }

    /**
     * @brief This function is called if the FMM is not computed on
     * all the standards levels
     *
     */
    void internal_M2L(){
570 571 572 573 574 575 576 577 578 579
        if(this->kernel->getKernelFct().m2l_ext){
            if(upperLimit > 1){ // if upperLimit == 1, then, M2L has been
                // done at level 2, and hence all the far
                // field has been calculated.
                //Starting at the lower level where the M2L has not been done.
                typename OctreeClass::Iterator octreeIterator(octree); //lvl : 1

                while(octreeIterator.level() != upperLimit){
                    octreeIterator.moveDown();
                }
580

581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627
                //I'm at the upperLimit, so the lowest level where M2L has been done.
                do{
                    CoreCell * currentTgt = octreeIterator.getCurrentCell(); // This one is targeted

                    //Then, we get the interaction list at this lvl. This will provide us with lots of source cells.
                    const CoreCell * currentInteractionList[343];
                    //Get an iterator for the sources
                    typename OctreeClass::Iterator upAndDownIterator = octreeIterator;

                    {//This is supposed to be done for multiple level. You
                        //need to go up until level 2. And then, to go down
                        //until level upperLimit. I think it's possible ...
                        while(upAndDownIterator.level() >= 2){
                            upAndDownIterator.moveUp();

                            //There, we get the interaction list of all parents of tgt cell
                            const int nbInteract = octree->getInteractionNeighbors(currentInteractionList,
                                                                                   upAndDownIterator.getCurrentGlobalCoordinate(),
                                                                                   upAndDownIterator.level());
                            int currentLevel = upAndDownIterator.level();
                            if(nbInteract){
                                //Then, we do M2L for each child at level upperLimit of each 343 Interaction cells.
                                for(int idxSrc = 0; idxSrc < 343 ; ++idxSrc){
                                    if(currentInteractionList[idxSrc]){//Check if it exist
                                        const CoreCell * currentSource = currentInteractionList[idxSrc]; //For clarity, will be otpimised out, anyway
                                        MortonIndex idx = currentSource->getMortonIndex();

                                        //At this point, we instanciate
                                        //the number of child needed.
                                        //This only depends on diffenrence
                                        //between current level and
                                        //upperLimit level
                                        int totalNumberOfChild = FMath::pow(8,upperLimit-currentLevel);

                                        for(int idxChildSrc = 0; idxChildSrc < totalNumberOfChild ; ++idxChildSrc){//For all 8^{number of levels to down} children
                                            MortonIndex indexOfChild = ((idx << 3*(upperLimit-currentLevel))+idxChildSrc);
                                            CoreCell * src = octree->getCell(indexOfChild,upperLimit); //Get the cell
                                            if(src){//check if it exists
                                                FTreeCoordinate srcCoord = src->getCoordinate();
                                                FTreeCoordinate tgtCoord = currentTgt->getCoordinate();
                                                //Build tree coord translation vector
                                                FTreeCoordinate transfer;
                                                transfer.setPosition(tgtCoord.getX()-srcCoord.getX(),
                                                                     tgtCoord.getY()-srcCoord.getY(),
                                                                     tgtCoord.getZ()-srcCoord.getZ());
                                                kernel->M2L_Extended(src,currentTgt,transfer,octreeIterator.level());
                                            }
628 629 630 631 632 633
                                        }
                                    }
                                }
                            }
                        }
                    }
634 635 636 637 638
                }while(octreeIterator.moveRight());
            }
            else{
                FAssertLF("No reasons to be there, seriously ...\nExiting anyway...");
            }
639
        }
640

641 642
    }

643
    virtual void execute_fmm(){
644
        FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
645
        switch(FScalFMMEngine<FReal>::Algorithm){
646 647
        case 0:
            {
648
                typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassSeq;
649 650
                AlgoClassSeq * algoSeq = new AlgoClassSeq(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoSeq;
651
                FScalFMMEngine<FReal>::abstrct = algoSeq;
652
                //algoSeq->execute(); will be done later
653 654 655 656 657
                break;
            }
        case 1:
            {
                typedef FFmmAlgorithmThread<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassThread;
658 659
                AlgoClassThread*  algoThread = new AlgoClassThread(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoThread;
660
                FScalFMMEngine<FReal>::abstrct = algoThread;
661
                //algoThread->execute(); will be done later
662 663 664 665
                break;
            }
        case 2:
            {
666
                typedef FFmmAlgorithmPeriodic<FReal,OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassPeriodic;
667 668 669 670 671
                AlgoClassPeriodic algoPeriod(octree,2);
                algoPeriod.setKernel(kernel);
                algoPeriod.execute();
                break;
            }
672 673
        case 3:
            {
674 675 676
                typedef FFmmAlgorithmThreadTsm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassTargetSource;
                AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoTS;
677
                FScalFMMEngine<FReal>::abstrct = algoTS;
678 679
                //algoTS->execute(); will be done later
                break;
680
            }
681
        default :
682
            std::cout<< "No algorithm found (probably for strange reasons) : "<< FScalFMMEngine<FReal>::Algorithm <<" exiting" << std::endl;
683
        }
684

685 686
        if (FScalFMMEngine<FReal>::Algorithm != 2){
            if(upperLimit != 2){
687
                (FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2M | FFmmM2M | FFmmM2L, upperLimit, treeHeight);
688 689 690
                printf("\tUpPass finished\n");
                internal_M2L();
                printf("\tStrange M2L finished\n");
691
                (FScalFMMEngine<FReal>::abstrct)->execute(FFmmL2L | FFmmL2P | FFmmP2P, upperLimit, treeHeight);
692 693 694
                printf("\tDownPass finished\n");
            }
            else{
695 696 697
                if(octree->getHeight() == 2){
                    (FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2P);
                }else{
698
                    (FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2M | FFmmM2M | FFmmM2L | FFmmL2L | FFmmL2P | FFmmP2P);
699
                }
700 701
            }
        }
702
    }
703

704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766
        virtual void execute_fmm_far_field(){
        FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
        switch(FScalFMMEngine<FReal>::Algorithm){
        case 0:
            {
                typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassSeq;
                AlgoClassSeq * algoSeq = new AlgoClassSeq(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoSeq;
                FScalFMMEngine<FReal>::abstrct = algoSeq;
                //algoSeq->execute(); will be done later
                break;
            }
        case 1:
            {
                typedef FFmmAlgorithmThread<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassThread;
                AlgoClassThread*  algoThread = new AlgoClassThread(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoThread;
                FScalFMMEngine<FReal>::abstrct = algoThread;
                //algoThread->execute(); will be done later
                break;
            }
        case 2:
            {
                typedef FFmmAlgorithmPeriodic<FReal,OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassPeriodic;
                AlgoClassPeriodic algoPeriod(octree,2);
                algoPeriod.setKernel(kernel);
                algoPeriod.execute();
                break;
            }
        case 3:
            {
                typedef FFmmAlgorithmThreadTsm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassTargetSource;
                AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
                FScalFMMEngine<FReal>::algoTimer = algoTS;
                FScalFMMEngine<FReal>::abstrct = algoTS;
                //algoTS->execute(); will be done later
                break;
            }
        default :
            std::cout<< "No algorithm found (probably for strange reasons) : "<< FScalFMMEngine<FReal>::Algorithm <<" exiting" << std::endl;
        }

        if (FScalFMMEngine<FReal>::Algorithm != 2){
            if(upperLimit != 2){
                (FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2M | FFmmM2M | FFmmM2L, upperLimit, treeHeight);
                printf("\tUpPass finished\n");
                internal_M2L();
                printf("\tStrange M2L finished\n");
                (FScalFMMEngine<FReal>::abstrct)->execute(FFmmL2L | FFmmL2P, upperLimit, treeHeight);
                printf("\tDownPass finished\n");
            }
            else{
                if(octree->getHeight() == 2){
                    //(FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2P);
                    std::cout<<"Nothing to be done\n";
                }else{
                    (FScalFMMEngine<FReal>::abstrct)->execute(FFmmP2M | FFmmM2M | FFmmM2L | FFmmL2L | FFmmL2P);
                }
            }
        }
    }


767
    virtual void intern_dealloc_handle(Callback_free_cell userDeallocator){
768
        free_cell(userDeallocator, ContainerClass::GetFreeLeaf());
769
    }
770 771 772 773
};


#endif //FUSERKERNELENGINE_HPP