FInterEngine.hpp 26 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
// ===================================================================================
// Copyright ScalFmm 2014 I
// 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".
// ===================================================================================


/**
 * @file This file contains a class that inherits from FScalFMMEngine,
 * and will implement the API functions for Interpolations kernels.
 */

#ifndef FINTERENGINE_HPP
#define FINTERENGINE_HPP

#include "FScalFMMEngine.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
27
//#include "Kernels/P2P/FP2PLeafInterface.hpp"
28 29
#include "Arranger/FOctreeArranger.hpp"
#include "Arranger/FArrangerPeriodic.hpp"
30
#include "Arranger/FBasicParticleContainerIndexedMover.hpp"
31

32 33
#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithm.hpp"
34
#include "Core/FFmmAlgorithmPeriodic.hpp"
35 36 37 38 39 40


/**
 * @class FInterEngine implements API for Interpolations kernels, its
 * templates can be ChebCell/ChebKernel or UnifCell/UnifKernel
 */
41
template<class FReal, class InterCell,class InterKernel,
42 43
         class ContainerClass = FP2PParticleContainerIndexed<FReal>,
         class LeafClass = FSimpleLeaf<FReal, FP2PParticleContainerIndexed<FReal> >,
44
         class MatrixKernelClass = FInterpMatrixKernelR<FReal> >
45 46
class FInterEngine : public FScalFMMEngine{
private:
47
    //Typedef on the octree class, in order to clarify following code
48
    typedef FOctree<FReal,InterCell,ContainerClass,LeafClass> OctreeClass;
49
    //typedef FP2PLeafInterface<OctreeClass>            LeafInterface;
50 51 52


    //Typedef on Octree Arranger, in order to clarify following code
53
    typedef FBasicParticleContainerIndexedMover<FReal,OctreeClass, ContainerClass> MoverClass;
54
    typedef FOctreeArranger<FReal,OctreeClass, ContainerClass, MoverClass> ArrangerClass;
55
    typedef FArrangerPeriodic<FReal,OctreeClass, ContainerClass, MoverClass> ArrangerClassPeriodic;
56

57 58
    //Pointer to the kernel to be executed
    InterKernel * kernel;
59
    MatrixKernelClass * matrix;
60
    //Link to the tree
61 62 63
    OctreeClass * octree;
    ArrangerClass * arranger;

64 65
public:
    /**
66 67
     * @brief Constructor : build the tree and the interpolation
     * kernel
68 69 70 71 72
     * @param TreeHeight Height of the tree
     * @param BoxWidth box Width
     * @param BoxCenter double[3] coordinate of the center of the
     * simulation box
     */
73
    FInterEngine(scalfmm_kernel_type KernelType) :
74
        kernel(nullptr), matrix(nullptr), octree(nullptr),arranger(nullptr){
75 76 77 78
        kernelType = KernelType;
    }

    void build_tree(int TreeHeight, double BoxWidth , double * BoxCenter,User_Scalfmm_Cell_Descriptor notUsedHere){
79
        octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint<FReal>(BoxCenter));
80
        this->matrix = new MatrixKernelClass();
81
        this->kernel = new InterKernel(TreeHeight,BoxWidth,FPoint<FReal>(BoxCenter),matrix);
82 83
    }

84

85 86
    //TODO free kernel too
    ~FInterEngine(){
87 88 89 90 91 92
        delete matrix;
        delete octree;
        delete kernel;
        if(arranger){
            delete arranger;
        }
93 94 95 96 97
    }

    //Inserting array of position
    void tree_insert_particles_xyz(int NbPositions, double * XYZ){
        for(int idPart = 0; idPart<NbPositions ; ++idPart){
98
            octree->insert(FPoint<FReal>(&XYZ[3*idPart]),idPart);
99 100 101 102 103 104 105
        }
        nbPart += NbPositions;
    }

    //Inserting arrayS of position
    void tree_insert_particles(int NbPositions, double * X, double * Y, double * Z){
        for(int idPart = 0; idPart<NbPositions ; ++idPart){
106
            octree->insert(FPoint<FReal>(X[idPart],Y[idPart],Z[idPart]),idPart);
107 108 109 110 111 112 113 114
        }
        nbPart += NbPositions;
    }

    //Set the physical values
    void set_physical_values(int nbPhysicalValues,double * physicalValues){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
115
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
116
                FSize nbPartThere = sources->getNbParticles();
117
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
118 119 120 121 122 123 124 125 126 127 128
                    sources->getPhysicalValues()[idxPart] = physicalValues[indexes[idxPart]];
                }
            });
    }

    //Set only a subpart of physical values
    //Algorithm : loop over each leaf, and then search in user array
    //if any index matches
    void set_physical_values_npart( int nbPhysicalValues, int* idxOfParticles, double * physicalValues){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
129
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
130
                FSize nbPartThere = sources->getNbParticles();
131
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbPhysicalValues && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            sources->getPhysicalValues()[idxPart] = physicalValues[indexes[idxPart]];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    //get back the physical values
    void get_physical_values( int nbPhysicalValues, double * physicalValues){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
151
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
152
                FSize nbPartThere = sources->getNbParticles();
153
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
154 155 156 157 158 159 160 161 162
                    physicalValues[indexes[idxPart]] = sources->getPhysicalValues()[idxPart];
                }
            });
    }

    //Same algorithm as in set_physical_values_npart
    void get_physical_values_npart( int nbPhysicalValues, int* idxOfParticles, double * physicalValues){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
163
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
164
                FSize nbPartThere = sources->getNbParticles();
165
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbPhysicalValues && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            physicalValues[indexes[idxPart]] = sources->getPhysicalValues()[idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    void get_forces_xyz( int nbParts, double * forcesToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
184
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
185
                FSize nbPartThere = sources->getNbParticles();
186
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
187 188 189 190 191 192 193 194 195 196
                    forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
                    forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
                    forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
                }
            });
    }

    void get_forces_xyz_npart(int nbParts, int* idxOfParticles , double * forcesToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
197
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
198
                FSize nbPartThere = sources->getNbParticles();
199
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbParts && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
                            forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
                            forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    void get_forces( int nbParts, double * fX, double* fY, double* fZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
220
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
221
                FSize nbPartThere = sources->getNbParticles();
222
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
223 224 225 226 227
                    fX[indexes[idxPart]] = sources->getForcesX()[idxPart];
                    fY[indexes[idxPart]] = sources->getForcesY()[idxPart];
                    fZ[indexes[idxPart]] = sources->getForcesZ()[idxPart];
                }
            });
228
    }
229 230 231 232

    void get_forces_npart(int nbParts, int* idxOfParticles ,double * fX, double* fY, double* fZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
233
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
234
                FSize nbPartThere = sources->getNbParticles();
235
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbParts && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            fX[indexes[idxPart]] = sources->getForcesX()[idxPart];
                            fY[indexes[idxPart]] = sources->getForcesY()[idxPart];
                            fZ[indexes[idxPart]] = sources->getForcesZ()[idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    //To set initial condition
254
    void set_forces_xyz( int nbParts, double * forcesToRead){
255 256
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
257
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
258
                FSize nbPartThere = sources->getNbParticles();
259
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
260 261 262
                    sources->getForcesX()[idxPart] = forcesToRead[indexes[idxPart]*3+0];
                    sources->getForcesY()[idxPart] = forcesToRead[indexes[idxPart]*3+1];
                    sources->getForcesZ()[idxPart] = forcesToRead[indexes[idxPart]*3+2];
263 264 265
                }
            });
    }
266
    void set_forces_xyz_npart( int nbParts, int* idxOfParticles, double * forcesToRead){
267 268
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
269
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
270
                FSize nbPartThere = sources->getNbParticles();
271
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
272 273 274 275
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbParts && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
276 277 278
                            sources->getForcesX()[idxPart] = forcesToRead[indexes[idxPart]*3+0];
                            sources->getForcesY()[idxPart] = forcesToRead[indexes[idxPart]*3+1];
                            sources->getForcesZ()[idxPart] = forcesToRead[indexes[idxPart]*3+2];
279 280 281 282 283 284 285 286 287 288 289 290
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }
    void set_forces( int nbParts, double * fX, double* fY, double* fZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
291
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
292
                FSize nbPartThere = sources->getNbParticles();
293
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
294 295 296 297 298 299 300 301 302
                    sources->getForcesX()[idxPart] = fX[indexes[idxPart]];
                    sources->getForcesY()[idxPart] = fY[indexes[idxPart]];
                    sources->getForcesZ()[idxPart] = fZ[indexes[idxPart]];
                }
            });
    }
    void set_forces_npart( int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
303
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
304
                FSize nbPartThere = sources->getNbParticles();
305
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbParts && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            sources->getForcesX()[idxPart] = fX[indexes[idxPart]];
                            sources->getForcesY()[idxPart] = fY[indexes[idxPart]];
                            sources->getForcesZ()[idxPart] = fZ[indexes[idxPart]];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    //Set the potentials
    void set_potentials(int nbPotentials,double * potentialsToRead){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
327
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
328
                FSize nbPartThere = sources->getNbParticles();
329
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
330 331 332 333 334 335 336 337 338 339 340
                    sources->getPotentials()[idxPart] = potentialsToRead[indexes[idxPart]];
                }
            });
    }

    //Set only a subpart of potentials
    //Algorithm : loop over each leaf, and then search in user array
    //if any index matches
    void set_potentials_npart( int nbPotentials, int* idxOfParticles, double * potentialsToRead){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
341
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
342
                FSize nbPartThere = sources->getNbParticles();
343
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbPotentials && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            sources->getPotentials()[idxPart] = potentialsToRead[indexes[idxPart]];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    //get back the potentials
    void get_potentials( int nbPotentials, double * potentialsToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
362
                ContainerClass * sources = leaf->getTargets();
363
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
364
                FSize nbPartThere = sources->getNbParticles();
365
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
366 367 368 369 370 371 372 373 374
                    potentialsToFill[indexes[idxPart]] = sources->getPotentials()[idxPart];
                }
            });
    }

    //Same algorithm as in set_potentials_npart
    void get_potentials_npart( int nbPotentials, int* idxOfParticles, double * potentialsToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
375
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
376
                FSize nbPartThere = sources->getNbParticles();
377
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < nbPotentials && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            potentialsToFill[indexes[idxPart]] = sources->getPotentials()[idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

393 394 395
    void get_positions_xyz(int NbPositions, double * positionsToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
396
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
397
                FSize nbPartThere = sources->getNbParticles();
398
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
399 400 401 402 403 404 405 406 407 408
                    positionsToFill[indexes[idxPart]*3+0] = sources->getPositions()[0][idxPart];
                    positionsToFill[indexes[idxPart]*3+1] = sources->getPositions()[1][idxPart];
                    positionsToFill[indexes[idxPart]*3+2] = sources->getPositions()[2][idxPart];
                }
            });
    }

    void get_positions_xyz_npart(int NbPositions, int * idxOfParticles, double * positionsToFill){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
409
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
410
                FSize nbPartThere = sources->getNbParticles();
411
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < NbPositions && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            positionsToFill[indexes[idxPart]*3+0] =  sources->getPositions()[0][idxPart];
                            positionsToFill[indexes[idxPart]*3+1] =  sources->getPositions()[1][idxPart];
                            positionsToFill[indexes[idxPart]*3+2] =  sources->getPositions()[2][idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }

    void get_positions( int NbPositions, double * X, double * Y , double * Z){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
432
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
433
                FSize nbPartThere = sources->getNbParticles();
434
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
435 436 437 438 439 440 441 442 443 444
                    X[indexes[idxPart]] = sources->getPositions()[0][idxPart];
                    Y[indexes[idxPart]] = sources->getPositions()[1][idxPart];
                    Z[indexes[idxPart]] = sources->getPositions()[2][idxPart];
                }
            });
    }

    void get_positions_npart(int NbPositions, int * idxOfParticles,double * X, double * Y , double * Z){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
445
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
446
                FSize nbPartThere = sources->getNbParticles();
447
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470
                    int iterPart = 0;
                    bool notFoundYet = true;
                    while(iterPart < NbPositions && notFoundYet){
                        if(indexes[idxPart] == idxOfParticles[iterPart]){
                            X[indexes[idxPart]] =  sources->getPositions()[0][idxPart];
                            Y[indexes[idxPart]] =  sources->getPositions()[1][idxPart];
                            Z[indexes[idxPart]] =  sources->getPositions()[2][idxPart];
                            notFoundYet = false;
                        }
                        else{
                            ++iterPart;
                        }
                    }
                }
            });
    }


    //Arranger parts : following function provide a way to move parts
    //inside the tree
    void add_to_positions_xyz(int NbPositions,double * updatedXYZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
471
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
472
                FSize nbPartThere = sources->getNbParticles();
473
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
474 475 476 477 478 479 480 481 482 483
                    sources->getWPositions()[0][idxPart] += updatedXYZ[indexes[idxPart]*3+0];
                    sources->getWPositions()[1][idxPart] += updatedXYZ[indexes[idxPart]*3+1];
                    sources->getWPositions()[2][idxPart] += updatedXYZ[indexes[idxPart]*3+2];
                }
            });
    }

    void add_to_positions(int NbPositions,double * X, double * Y , double * Z){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
484
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
485
                FSize nbPartThere = sources->getNbParticles();
486
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
487 488 489 490 491 492 493 494 495 496 497
                    sources->getWPositions()[0][idxPart] += X[indexes[idxPart]];
                    sources->getWPositions()[1][idxPart] += Y[indexes[idxPart]];
                    sources->getWPositions()[2][idxPart] += Z[indexes[idxPart]];
                }
            });
    }


    void set_positions_xyz(int NbPositions, double * updatedXYZ){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
498
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
499
                FSize nbPartThere = sources->getNbParticles();
500
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
501 502 503 504 505 506 507 508 509 510
                    sources->getWPositions()[0][idxPart] = updatedXYZ[indexes[idxPart]*3+0];
                    sources->getWPositions()[1][idxPart] = updatedXYZ[indexes[idxPart]*3+1];
                    sources->getWPositions()[2][idxPart] = updatedXYZ[indexes[idxPart]*3+2];
                }
            });
    }

    void set_positions(int NbPositions, double * X, double * Y, double * Z){
        octree->forEachLeaf([&](LeafClass* leaf){
                ContainerClass * sources = leaf->getSrc();
511
                const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
512
                FSize nbPartThere = sources->getNbParticles();
513
                for(FSize idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
514 515 516 517 518 519 520
                    sources->getWPositions()[0][idxPart] = X[indexes[idxPart]];
                    sources->getWPositions()[1][idxPart] = Y[indexes[idxPart]];
                    sources->getWPositions()[2][idxPart] = Z[indexes[idxPart]];
                }
            });
    }

521
    //Simple call to FScalFMMEngine method with good template
522
    void reset_tree(Callback_reset_cell /*not used*/){
523 524 525
        generic_reset_tree<FReal,ContainerClass,InterCell,LeafClass>(octree);
    }

526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
    void update_tree(){
        if(arranger){
            arranger->rearrange();
        }
        else{
            if(Algorithm == 2){ //case in wich the periodic algorithm is used
                arranger = new ArrangerClassPeriodic(octree);
                arranger->rearrange();
            }
            else{
                arranger = new ArrangerClass(octree);
                arranger->rearrange();
            }
        }
    }
541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558


    void execute_fmm(){
        switch(Algorithm){
        case 0:
            {
                typedef FFmmAlgorithm<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassSeq;
                AlgoClassSeq algoSeq(octree,kernel);
                algoSeq.execute();
                break;
            }
        case 1:
            {
                typedef FFmmAlgorithmThread<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassThread;
                AlgoClassThread algoThread(octree,kernel);
                algoThread.execute();
                break;
            }
559 560
        case 2:
            {
561
                typedef FFmmAlgorithmPeriodic<FReal,OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassPeriodic;
562
                AlgoClassPeriodic algoPeriod(octree,2);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
563
                algoPeriod.setKernel(kernel);
564 565 566
                algoPeriod.execute();
                break;
            }
567 568 569 570 571
        default :
            std::cout<< "No algorithm found (probably for strange reasons) : "<< Algorithm <<" exiting" << std::endl;
        }
    }

572 573 574
    void intern_dealloc_handle(Callback_free_cell unUsed){
        //this->~FInterEngine();
    }
575 576 577 578
};


#endif