FScalfmmApiInit.cpp 31.7 KB
Newer Older
1 2 3 4 5 6
/** It should be compiled with C export */
extern "C" {
#include "CScalfmmApi.h"
}

#include "FInterEngine.hpp"
7
#include "FUserKernelEngine.hpp"
8

9 10 11 12 13 14 15 16 17

/**
 * Define here static member
 */
Scalfmm_Cell_Descriptor CoreCell::user_cell_descriptor;
template<class FReal>
Scalfmm_Leaf_Descriptor FUserLeafContainer<FReal>::user_leaf_descriptor;


18 19 20
extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,
                                         double* BoxCenter, */
                                       scalfmm_kernel_type KernelType,
21 22
                                       scalfmm_algorithm algo){
    ScalFmmCoreHandle<double> * handle = new ScalFmmCoreHandle<double>();
23
    typedef double FReal;
24

25
    if(algo == source_target){
26

27 28
        switch(KernelType){
        case 0:
29 30
            typedef FUserLeafContainer<FReal>                                               UserContainerClass;
            typedef FTypedLeaf<FReal,UserContainerClass>                                         UserLeafClass;
31

32
            handle->engine = new FUserKernelEngine<FReal,UserLeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType,algo);
33
            break;
34

35 36 37 38 39
        case 1:
            //TODO typedefs
            typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
            typedef FTypedChebCell<FReal,7>                                                   ChebCell;
            typedef FTypedLeaf<FReal,ContainerClass>                                         LeafClass;
40

41 42
            typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
            typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
43

44
            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType,algo);
45 46 47 48 49
            break;
            // case 2:
            //     //TODO typedefs
            //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
            //     typedef FUnifCell<7>                                                         UnifCell;
50

51 52 53 54 55 56 57 58 59 60 61 62
            //     typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
            //     typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7>           UnifKernel;

            //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
            //     break;

        default:
            std::cout<< "Kernel type unsupported" << std::endl;
            exit(0);
            break;
        }
    }
63 64 65 66 67 68 69 70
    else{
        // if(algo == adaptiv){
        //     //Temporary
        //     handle->engine = new FAdaptEngine<FReal,4>(KernelType,algo);
        //}else{
            //No Source/Targets distinction
            switch(KernelType){
            case 0:
71 72
                typedef FUserLeafContainer<FReal>                            UserContainerClass;
                typedef FSimpleLeaf<FReal,UserContainerClass>                         UserLeafClass;
73

74
                handle->engine = new FUserKernelEngine<FReal,UserLeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType,algo);
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
                break;

            case 1:
                //TODO typedefs
                typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
                //typedef FChebCell<FReal,7>                                                   ChebCell;
                typedef FTypedChebCell<FReal,7>                                                   ChebCell;
                typedef FSimpleLeaf<FReal,ContainerClass>                                         LeafClass;

                typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
                typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;

                handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType,algo);
                break;
                // case 2:
                //     //TODO typedefs
                //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
                //     typedef FUnifCell<7>                                                         UnifCell;

                //     typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
                //     typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7>           UnifKernel;

                //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
                //     break;
            default:
                std::cout<< "Kernel type unsupported" << std::endl;
                exit(0);
                break;
            }
            // }
105
    }
106
   return handle;
107 108
}

109
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
110 111 112
    ((ScalFmmCoreHandle<double> *) handle)->engine->intern_dealloc_handle(userDeallocator);
    delete ((ScalFmmCoreHandle<double> *) handle)->engine ;
    delete (ScalFmmCoreHandle<double> *) handle;
113
}
114

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

/**
 * @brief Init function for distributed version (MPI).
 *
 */
#ifdef SCALFMM_USE_MPI

#include "Utils/FMpi.hpp"
#include "FUserKernelDistrEngine.hpp"

extern "C" scalfmm_handle scalfmm_init_distributed( scalfmm_kernel_type KernelType,
                                                    scalfmm_algorithm algo,
                                                    const MPI_Comm comm){

    ScalFmmCoreHandle<double> * handle = new ScalFmmCoreHandle<double>();

    //Only the User Defined Kernel version (UDK) is available.
    typedef double FReal;
133 134
    typedef FUserLeafContainer<FReal>                            ContainerClass;
    typedef FSimpleLeaf<FReal,ContainerClass>                         LeafClass;
135 136 137 138 139 140 141

    handle->engine = (new FUserKernelDistrEngine<FReal,LeafClass>(KernelType,algo,comm));
    return handle;
}

#endif

142 143 144 145 146 147
/**
 * This parts implements all the function defined in FChebInterface.h
 * using the Chebyshev classes
 */
#ifndef CHEBINTERFACE_HPP
#define CHEBINTERFACE_HPP
148 149 150 151

#warning "Compiling Cheb Interface"


152 153 154 155 156 157 158 159 160 161 162
extern "C" {
#include "Kernels/Chebyshev/FChebInterface.h"
}


#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Components/FSimpleLeaf.hpp"


typedef struct FChebCell_struct{
    //Store what's needed
163
    FChebCell<double,7> * cell;
164 165
}ChebCellStruct;

166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
typedef struct FChebLeaf_struct{
    //Store a P2PParticleContainer
    FP2PParticleContainerIndexed<double>* container;
}ChebLeafStruct;


//Initialize leaves
extern "C" ChebLeafStruct * ChebLeafStruct_create(FSize nbPart){
    FP2PParticleContainerIndexed<double>* newCont = new FP2PParticleContainerIndexed<double>();
    newCont->reserve(nbPart);
    ChebLeafStruct * newLeaf = new ChebLeafStruct();
    newLeaf->container = newCont;
    return newLeaf;
}

//Delete leaves
extern "C" void ChebLeafStruct_free(void * leafData){
    ChebLeafStruct * leaf = reinterpret_cast<ChebLeafStruct *>(leafData);
    delete leaf->container;
    delete leaf;
}

//Fill leaves once partitionning is done
extern "C" void ChebLeafStruct_fill(FSize nbPart, const FSize * idxPart,
                                    long long morton_index, void * leafData,
                                    void * userData){
    UserData * userDataKernel = reinterpret_cast<UserData *>(userData);
    ChebLeafStruct * leaf = reinterpret_cast<ChebLeafStruct *>(leafData);
    FP2PParticleContainerIndexed<double>* container = leaf->container;

    for(int i=0 ; i<nbPart; ++i){
        FPoint<double> pos{userDataKernel->insertedPositions[idxPart[i]*3 + 0],
                userDataKernel->insertedPositions[idxPart[i]*3 + 1],
                userDataKernel->insertedPositions[idxPart[i]*3 + 2]};
        double phi = userDataKernel->myPhyValues[idxPart[i]];
        container->push(pos,idxPart[i],phi);
    }
}


extern "C" void ChebLeafStruct_get_back_results(void * leafData,
                                                double ** forceXptr,  double ** forceYptr,  double ** forceZptr,
                                                double ** potentialsptr){
    ChebLeafStruct * leaf = reinterpret_cast<ChebLeafStruct *>(leafData);
    *forceXptr = leaf->container->getForcesX();
    *forceYptr = leaf->container->getForcesY();
    *forceZptr = leaf->container->getForcesZ();
    *potentialsptr = leaf->container->getPotentials();
}
215 216 217 218

//How to create/destroy cells
extern "C" ChebCellStruct * ChebCellStruct_create(long long int inIndex,int * position){
    ChebCellStruct * newCell = new ChebCellStruct();
219
    newCell->cell = new FChebCell<double,7>();
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235
    newCell->cell->setMortonIndex(inIndex);
    newCell->cell->setCoordinate(position[0],position[1],position[2]);
    return newCell;
}

extern "C" void ChebCellStruct_free(ChebCellStruct * inCell){
    if(inCell->cell) {
        delete inCell->cell;
    }
    delete inCell;
}



typedef struct FChebKernel_struct{
    //Be ready full duplication go there !!!
236 237
    FChebSymKernel<double,FChebCell<double,7>,FP2PParticleContainerIndexed<double>,FInterpMatrixKernelR<double>,7> ** kernel;
    FInterpMatrixKernelR<double> * matrix;
238 239
} ChebKernelStruct;

240 241


242 243 244 245
//Kernel functions
extern "C" ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
                                                      double inBoxWidth,
                                                      double* inBoxCenter){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
246
    //typedef to lighten the writing
247
    typedef FChebSymKernel<double,FChebCell<double,7>, FP2PParticleContainerIndexed<double>, FInterpMatrixKernelR<double>,7> KernelClass;
248
    ChebKernelStruct * newKernel = new ChebKernelStruct();
249
    newKernel->matrix= new FInterpMatrixKernelR<double>();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
250 251 252 253
    int nb_threads = omp_get_max_threads();
    newKernel->kernel = new KernelClass*[nb_threads];
    newKernel->kernel[0]=new KernelClass(inTreeHeight,
                                         inBoxWidth,
254
                                         FPoint<double>(inBoxCenter[0], inBoxCenter[1], inBoxCenter[2]),
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
255 256 257 258 259
                                         newKernel->matrix);

    for(int idThreads=1 ; idThreads<nb_threads ; ++idThreads){
        newKernel->kernel[idThreads] = new KernelClass(*newKernel->kernel[0]);
    }
260 261 262 263 264
    return newKernel;
}

extern "C" void ChebKernelStruct_free(void *inKernel){
    delete reinterpret_cast<ChebKernelStruct *>(inKernel)->matrix;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
265 266 267 268 269
    int nb_threads = omp_get_max_threads();
    for(int idThreads=0 ; idThreads<nb_threads ; ++idThreads){
        delete reinterpret_cast<ChebKernelStruct *>(inKernel)->kernel[idThreads];
    }
    delete [] reinterpret_cast<ChebKernelStruct *>(inKernel)->kernel;
270 271 272
    delete reinterpret_cast<ChebKernelStruct *>(inKernel);
}

273 274 275 276 277 278 279
/**
* New one, intended to work with internal FP2PParticleContainer instead of filling a new one
*/
extern "C" void ChebKernel_P2M(void * leafCell, void * leafData, FSize nbParticles,
                               const FSize *particleIndexes, void * inKernel){
    //get the real leaf
    ChebLeafStruct * leaf = reinterpret_cast<ChebLeafStruct *>(leafData);
280

281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
    //get the real cell struct
    ChebCellStruct * realCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    FChebCell<double,7> * realCell = realCellStruct->cell;

    //Identify thread number
    int id_thread = omp_get_thread_num();

    //get the real chebyshev struct
    UserData * userDataKernel = reinterpret_cast<UserData *>(inKernel);
    ChebKernelStruct * realKernel = userDataKernel->kernelStruct;

    realKernel->kernel[id_thread]->P2M(realCell, leaf->container);
}

extern "C" void ChebKernel_P2M_old(void * leafCell, void * leafData, FSize nbParticles,
                                   const FSize *particleIndexes, void * inKernel){
297
    //make temporary array of parts
298
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
299
    tempContainer->reserve(nbParticles);
300
    FPoint<double> pos;
301
    for(int i=0 ; i<nbParticles ; ++i){
302
        pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
303 304
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
305 306 307 308 309 310
        double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];

        tempContainer->push(pos,particleIndexes[i],Phi);
    }
    //get the real cell struct
    ChebCellStruct * realCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
311
    FChebCell<double,7> * realCell = realCellStruct->cell;
312

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
313 314 315
    //Identify thread number
    int id_thread = omp_get_thread_num();

316 317 318 319
    //get the real chebyshev struct
    UserData * userDataKernel = reinterpret_cast<UserData *>(inKernel);
    ChebKernelStruct * realKernel = userDataKernel->kernelStruct;

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
320
    realKernel->kernel[id_thread]->P2M(realCell, tempContainer);
321 322 323
    delete tempContainer;
}

324

325 326 327 328 329
extern "C" void  ChebKernel_M2M(int level, void* parentCell, int childPosition, void *childCell, void *inKernel){
    //Get our structures
    ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
    ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
    //get real cheb cell
330 331
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
332

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
333 334 335
    //Identify thread number
    int id_thread = omp_get_thread_num();

336 337
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
338
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyM2M(childPosition,
339 340
                                                                        childChebCell->getMultipole(0),
                                                                        parentChebCell->getMultipole(0));
341 342
}

343 344
extern "C" void ChebKernel_M2L(int level, void* targetCell, const int*neighborPositions,int size,
                               void** sourceCell, void* inKernel){
345 346 347
    //Get our structures
    ChebCellStruct * targetCellStruct = reinterpret_cast<ChebCellStruct *>(targetCell);
    //get real cheb cell
348
    FChebCell<double,7>* const targetChebCell = targetCellStruct->cell;
349 350

    //copy to an array of FChebCell
351 352 353
    FChebCell<double,7> const ** arrayOfChebCell = new const FChebCell<double,7>*[size];
    for(int i=0; i<size ; ++i){
        arrayOfChebCell[i] = reinterpret_cast<ChebCellStruct*>(sourceCell[i])->cell;
354
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
355 356 357 358

    //Identify thread number
    int id_thread = omp_get_thread_num();

359 360
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
361 362
    inKernelStruct->kernel[id_thread]->M2L(targetChebCell,arrayOfChebCell,neighborPositions,size,level);
    delete [] arrayOfChebCell;
363 364 365
}

extern "C" void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* inKernel){
366
    //Get our structures
367 368 369
    ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
    ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
    //get real cheb cell
370 371
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
372

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
373 374 375
    //Identify thread number
    int id_thread = omp_get_thread_num();

376 377
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
378
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyL2L(childPosition,
379 380
                                                                        parentChebCell->getLocal(0),
                                                                        childChebCell->getLocal(0));
381 382
}

383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
extern "C" void ChebKernel_L2P(void* leafCell, void * leafData, FSize nbParticles,
                               const FSize* particleIndexes, void* inKernel){
    //get the real leaf
    ChebLeafStruct * leaf = reinterpret_cast<ChebLeafStruct *>(leafData);

    //get the real cell struct
    ChebCellStruct * realCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    FChebCell<double,7> * realCell = realCellStruct->cell;

    //Identify thread number
    int id_thread = omp_get_thread_num();

    //get the real chebyshev struct
    UserData * userDataKernel = reinterpret_cast<UserData *>(inKernel);
    ChebKernelStruct * realKernel = userDataKernel->kernelStruct;

    realKernel->kernel[id_thread]->L2P(realCell,leaf->container);
}


extern "C" void ChebKernel_L2P_old(void* leafCell, void * leafData,FSize nbParticles,
                                   const FSize* particleIndexes, void* inKernel){
405
    //Create temporary FSimpleLeaf
406
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
407
    tempContainer->reserve(nbParticles);
408
    FPoint<double> pos;
409
    for(int i=0 ; i<nbParticles ; ++i){
410 411 412 413 414 415
        pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
        double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
        tempContainer->push(pos,particleIndexes[i],Phi);
    }
416 417 418
    //Get our structures
    ChebCellStruct * leafCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    //get real cheb cell
419
    FChebCell<double,7>* leafChebCell = leafCellStruct->cell;
420

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
421 422 423
    //Identify thread number
    int id_thread = omp_get_thread_num();

424 425 426
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
427
    inKernelStruct->kernel[id_thread]->L2P(leafChebCell,tempContainer);
428 429

    //Then retrieve the results
430 431 432
    double * forcesToFill     = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
    double * potentialsToFill = reinterpret_cast<UserData *>(inKernel)->potentials[id_thread];

433
    const FVector<FSize>& indexes = tempContainer->getIndexes();
434
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
435 436 437
        forcesToFill[indexes[idxPart]*3+0] += tempContainer->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContainer->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContainer->getForcesZ()[idxPart];
438
        potentialsToFill[indexes[idxPart]] += tempContainer->getPotentials()[idxPart];
439
    }
440 441 442 443 444

    delete tempContainer;
    tempContainer=nullptr;
}

445 446
void ChebKernel_P2P(void * targetleaf, FSize nbParticles, const FSize* particleIndexes, void ** sourceLeaves,
                    const FSize ** sourceParticleIndexes, FSize* sourceNbPart,
447
                    const int * sourcePosition,const int size, void* inKernel){
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481
    //First step, convert the target leaf in a P2P Particle container
    ChebLeafStruct * tgtLeaf = reinterpret_cast<ChebLeafStruct *>(targetleaf);

    //Same with the sources leaves
    std::vector<FP2PParticleContainerIndexed<double> *> arraySrcLeaves;
    arraySrcLeaves.reserve(size);
    for(int i=0 ; i<size ; ++i){
        if(sourceNbPart[i] != 0){
            arraySrcLeaves.push_back(reinterpret_cast<ChebLeafStruct *>(sourceLeaves[i])->container);
        }else{
            arraySrcLeaves.push_back(nullptr);
        }
    }

    //Then call P2P ?
    //Identify thread number
    int id_thread = omp_get_thread_num();

    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;

    //Empty tree coordinate
    int coord[3] = {0,0,0};

    inKernelStruct->kernel[id_thread]->P2P(FTreeCoordinate(coord),tgtLeaf->container,tgtLeaf->container,
                                           arraySrcLeaves.data(),sourcePosition,size);

}


void ChebKernel_P2P_old(void * targetLeaf, FSize nbParticles, const FSize* particleIndexes,
                        void ** sourceLeaves,
                        const FSize ** sourceParticleIndexes,FSize* sourceNbPart,
                        const int * sourcePosition,const int size, void* inKernel){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
482

483
    //Create temporary FSimpleLeaf for target
484
    FP2PParticleContainerIndexed<double>* tempContTarget = new FP2PParticleContainerIndexed<double>();
485 486
    tempContTarget->reserve(nbParticles);
    for(int i=0 ; i<nbParticles ; ++i){
487
        FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
488 489
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
490 491 492 493
        double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
        tempContTarget->push(pos,particleIndexes[i],Phi);
    }

494 495 496
    //Create size FSimpleLeaf for the size sources
    FP2PParticleContainerIndexed<double>** tempContSources = new FP2PParticleContainerIndexed<double>*[size];
    for(int idSource=0; idSource<size ; ++idSource){
497 498
        if(sourceNbPart[idSource] != 0){
            //Create container
499
            tempContSources[idSource] = new FP2PParticleContainerIndexed<double>();
500 501 502
            //Allocate memory
            tempContSources[idSource]->reserve(sourceNbPart[idSource]);
            //Store a ptr to the indices of that source leaf
503
            const FSize * indSource = sourceParticleIndexes[idSource];
504 505
            //Then, for each part in this source
            for(int i=0 ; i<sourceNbPart[idSource] ; ++i){
506
                FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3  ],
507 508
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+1],
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+2]);
509 510 511
                double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[indSource[i]];
                tempContSources[idSource]->push(pos,indSource[i],Phi);
            }
512
        } else{
513 514 515 516 517
            tempContSources[idSource] = nullptr;
        }
    }
    //Everything is fine, now, call Chebyshev P2P

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
518 519 520
    //Identify thread number
    int id_thread = omp_get_thread_num();

521 522 523 524 525 526
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;

    //Empty tree coordinate
    int coord[3] = {0,0,0};

527
    inKernelStruct->kernel[id_thread]->P2P(FTreeCoordinate(coord),tempContTarget,tempContTarget,tempContSources,sourcePosition,size);
528

529
    //get back forces & potentials
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
530
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
531 532
    double * potentialsToFill = reinterpret_cast<UserData *>(inKernel)->potentials[id_thread];

533
    const FVector<FSize>& indexes = tempContTarget->getIndexes();
534
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
535 536 537
        forcesToFill[indexes[idxPart]*3+0] += tempContTarget->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContTarget->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContTarget->getForcesZ()[idxPart];
538
        potentialsToFill[indexes[idxPart]] += tempContTarget->getPotentials()[idxPart];
539 540 541 542
    }

    //Note that sources are also modified.
    //get back sources forces
543
    for(int idSource = 0 ; idSource < size ; ++idSource){
544 545
        const FVector<FSize>& indexesSource = tempContSources[idSource]->getIndexes();
        const FSize nbPartInSource = sourceNbPart[idSource];
546
        for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
547 548 549
            forcesToFill[indexesSource[idxSourcePart]*3+0] += tempContSources[idSource]->getForcesX()[idxSourcePart];
            forcesToFill[indexesSource[idxSourcePart]*3+1] += tempContSources[idSource]->getForcesY()[idxSourcePart];
            forcesToFill[indexesSource[idxSourcePart]*3+2] += tempContSources[idSource]->getForcesZ()[idxSourcePart];
550
            potentialsToFill[indexesSource[idxSourcePart]] += tempContSources[idSource]->getPotentials()[idxSourcePart];
551 552 553 554
        }
    }

    //Release memory
555
    for(int idSource=0; idSource<size ; ++idSource){
556 557
        if(tempContSources[idSource]) delete tempContSources[idSource];
    }
558 559
    delete    tempContTarget;
    delete [] tempContSources;
560 561
}

562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596

void ChebKernel_P2PRemote(void * targetLeaf,FSize nbParticles, const FSize* particleIndexes,
                          void ** sourceLeaves,
                          const FSize ** sourceParticleIndexes,FSize * sourceNbPart,
                          const int * sourcePosition, const int size, void* inKernel){
    //First step, convert the target leaf in a P2P Particle container
    ChebLeafStruct * tgtLeaf = reinterpret_cast<ChebLeafStruct *>(targetLeaf);

    //Same with the sources leaves
    std::vector<FP2PParticleContainerIndexed<double> *> arraySrcLeaves;
    arraySrcLeaves.reserve(size);
    for(int i=0 ; i<size ; ++i){
        if(sourceNbPart[i] != 0){
            arraySrcLeaves.push_back(reinterpret_cast<ChebLeafStruct *>(sourceLeaves[i])->container);
        }else{
            arraySrcLeaves.push_back(nullptr);
        }
    }

    //Then call P2P ?
    //Identify thread number
    int id_thread = omp_get_thread_num();

    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;

    //Empty tree coordinate
    int coord[3] = {0,0,0};

    inKernelStruct->kernel[id_thread]->P2PRemote(FTreeCoordinate(coord),tgtLeaf->container,tgtLeaf->container,
                                                 arraySrcLeaves.data(),sourcePosition,size);

}


597
void ChebCell_reset(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell,void * inKernel){
598 599 600 601
    ChebCellStruct *  cellStruct = reinterpret_cast<ChebCellStruct *>(userCell);
    FChebCell<double,7>* chebCell = cellStruct->cell;
    chebCell->resetToInitialState();
}
602

603
FSize ChebCell_getSize(int level, long long morton_index){
604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676
    //Create fake one and ask for its size
    FChebCell<double,7>* chebCell = new FChebCell<double,7>();
    //then return what size is needed to store a cell
    FSize sizeToReturn = chebCell->getSavedSize();
    delete chebCell;
    return sizeToReturn;
}

/**
 * This is what the memory looks like
 * : [Morton]    [Coord][Multipole]        [Local]
 * : [1*longlong][3*int][double*VectorSize][double*VectorSize]
 */
void ChebCell_copy(void * userDatas, FSize size, void * memoryAllocated){
    //First cast inside outr struct
    ChebCellStruct *  cellStruct = reinterpret_cast<ChebCellStruct *>(userDatas);
    //Then get the FChebCell
    FChebCell<double,7>* chebCell = cellStruct->cell;
    //I know for sure there is enough place inside memory Allocated
    FSize cursor = 0;

    char * toWrite = static_cast<char * >(memoryAllocated);

    //Morton first
    long long here;
    here = chebCell->getMortonIndex();
    memcpy(&toWrite[cursor],&here, sizeof(long long));
    cursor += sizeof(long long);

    //FTreeCordinate then
    int coord[3] = {chebCell->getCoordinate().getX(),
                    chebCell->getCoordinate().getY(),
                    chebCell->getCoordinate().getZ()};
    memcpy(&toWrite[cursor],coord, sizeof(int)*3);
    cursor += 3*sizeof(int);
    //Upward datas
    memcpy(&toWrite[cursor],chebCell->getMultipole(0),chebCell->getVectorSize()*sizeof(double));
    cursor += sizeof(double)*chebCell->getVectorSize();
    //Downward datas
    memcpy(&toWrite[cursor],chebCell->getLocal(0),chebCell->getVectorSize()*sizeof(double));
    cursor += sizeof(double)*chebCell->getVectorSize();
}

/**
 * This is what the memory looks like
 * : [Morton]    [Coord][Multipole]        [Local]
 * : [1*longlong][3*int][double*VectorSize][double*VectorSize]
 */
void* ChebCell_restore(int level, void * arrayTobeRead){
    FSize cursor = 0;
    //First read Morton
    long long morton = (static_cast<long long *>(arrayTobeRead))[0];
    cursor += sizeof(long long);
    //Then Coord :
    int * coord = reinterpret_cast<int*>(&(static_cast<char*>(arrayTobeRead))[cursor]);
    cursor += sizeof(int)*3;

    //Create target Cell and Struct
    ChebCellStruct *  cellStruct = ChebCellStruct_create(morton,coord);

    //Then copy inside this cell the Multipole and the Local
    double * mult = reinterpret_cast<double*>(&(static_cast<char*>(arrayTobeRead))[cursor]);
    memcpy((cellStruct->cell)->getMultipole(0),mult,((cellStruct->cell)->getVectorSize())*sizeof(double));
    cursor += (cellStruct->cell)->getVectorSize()*sizeof(double);

    double * local = reinterpret_cast<double*>(&(static_cast<char*>(arrayTobeRead))[cursor]);
    memcpy((cellStruct->cell)->getLocal(0),local,((cellStruct->cell)->getVectorSize())*sizeof(double));
    cursor += (cellStruct->cell)->getVectorSize()*sizeof(double);

    //Yeah, can return, finally !!
    return cellStruct;
}

677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 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
/**
 * @brief Those fucntion implements the copy,restore and get_size
 * functions for leaves
 */
FSize ChebLeaf_getSize(FSize nbPart){
    //Test where we do not use leafData
    FP2PParticleContainerIndexed<double>* newCont = new FP2PParticleContainerIndexed<double>();
    newCont->reserve(nbPart);
    //fake push empty parts
    for(int i=0 ; i<nbPart ; ++i){
        FPoint<double> pos{0,0,0};
        newCont->push(pos,0,0,0,0,0);
    }

    FSize res = newCont->getSavedSize();
    delete newCont;

    return res;
}

/**
 * Copy Order : Positions XYZ, then attributes
 */
void ChebLeaf_copy(FSize nbPart,void * userLeaf, void * memAllocated){
    ChebLeafStruct * tgtLeaf = reinterpret_cast<ChebLeafStruct *>(userLeaf);
    FP2PParticleContainerIndexed<double>* container = tgtLeaf->container;

    char * toWrite = reinterpret_cast<char*>(memAllocated);

    FSize cursor = 0;
    //Loop over dimensions
    for(int i=0 ; i<3 ; ++i){
        memcpy(&toWrite[cursor],container->getPositions()[i],nbPart*sizeof(double));
        cursor += nbPart*sizeof(double);
    }

    //Then store attributes
    //First physical values
    memcpy(&toWrite[cursor],container->getPhysicalValuesArray(),sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fx
    memcpy(&toWrite[cursor],container->getForcesXArray(),sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fy
    memcpy(&toWrite[cursor],container->getForcesYArray(),sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fz
    memcpy(&toWrite[cursor],container->getForcesZArray(),sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;

    //Finally potentials
    memcpy(&toWrite[cursor],container->getPotentialsArray(),sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
}

void * ChebLeaf_restore(FSize nbPart,void * memToRead){
    ChebLeafStruct * tgtLeaf = ChebLeafStruct_create(nbPart);
    FP2PParticleContainerIndexed<double>* container = tgtLeaf->container;

    char * toRead = reinterpret_cast<char*>(memToRead);

    FSize cursor = 0;

    //Loop over dimensions
    for(int i=0 ; i<3 ; ++i){
        memcpy(container->getWPositions()[i],&toRead[cursor],nbPart*sizeof(double));
        cursor += nbPart*sizeof(double);
    }

    //Then store attributes
    //First physical values
    memcpy(container->getPhysicalValuesArray(),&toRead[cursor],sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fx
    memcpy(container->getForcesXArray(),&toRead[cursor],sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fy
    memcpy(container->getForcesYArray(),&toRead[cursor],sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;
    //Then fz
    memcpy(container->getForcesZArray(),&toRead[cursor],sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;

    //Finally potentials
    memcpy(container->getPotentialsArray(),&toRead[cursor],sizeof(double)*nbPart);
    cursor += sizeof(double)*nbPart;

    return tgtLeaf;
}

767
#endif