FScalfmmApiInit.cpp 18.3 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
extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double* BoxCenter, */scalfmm_kernel_type KernelType,
                                       scalfmm_algorithm algo){
    ScalFmmCoreHandle<double> * handle = new ScalFmmCoreHandle<double>();
13
    typedef double FReal;
14

15
    if(algo == source_target){
16

17 18 19 20
        switch(KernelType){
        case 0:
            typedef FP2PParticleContainerIndexed<FReal>           ContainerClass;
            typedef FTypedLeaf<FReal,ContainerClass>                                         LeafClass;
21

22 23
            handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
            break;
24

25 26 27 28 29
        case 1:
            //TODO typedefs
            typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
            typedef FTypedChebCell<FReal,7>                                                   ChebCell;
            typedef FTypedLeaf<FReal,ContainerClass>                                         LeafClass;
30

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

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

41 42 43 44 45 46 47 48 49 50 51 52
            //     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;
        }
    }
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
    else{
        // if(algo == adaptiv){
        //     //Temporary
        //     handle->engine = new FAdaptEngine<FReal,4>(KernelType,algo);
        //}else{
            //No Source/Targets distinction
            switch(KernelType){
            case 0:
                typedef FP2PParticleContainerIndexed<FReal>                            ContainerClass;
                typedef FSimpleLeaf<FReal,ContainerClass>                                   LeafClass;

                handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
                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;
            }
            // }
95
    }
96
   return handle;
97 98
}

99
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
100 101 102
    ((ScalFmmCoreHandle<double> *) handle)->engine->intern_dealloc_handle(userDeallocator);
    delete ((ScalFmmCoreHandle<double> *) handle)->engine ;
    delete (ScalFmmCoreHandle<double> *) handle;
103
}
104 105 106 107 108 109 110

/**
 * This parts implements all the function defined in FChebInterface.h
 * using the Chebyshev classes
 */
#ifndef CHEBINTERFACE_HPP
#define CHEBINTERFACE_HPP
111 112 113 114

#warning "Compiling Cheb Interface"


115 116 117 118 119 120 121 122 123 124 125
extern "C" {
#include "Kernels/Chebyshev/FChebInterface.h"
}


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


typedef struct FChebCell_struct{
    //Store what's needed
126
    FChebCell<double,7> * cell;
127 128 129 130 131 132
}ChebCellStruct;


//How to create/destroy cells
extern "C" ChebCellStruct * ChebCellStruct_create(long long int inIndex,int * position){
    ChebCellStruct * newCell = new ChebCellStruct();
133
    newCell->cell = new FChebCell<double,7>();
134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
    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 !!!
150 151
    FChebSymKernel<double,FChebCell<double,7>,FP2PParticleContainerIndexed<double>,FInterpMatrixKernelR<double>,7> ** kernel;
    FInterpMatrixKernelR<double> * matrix;
152 153 154 155 156 157
} ChebKernelStruct;

//Kernel functions
extern "C" ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
                                                      double inBoxWidth,
                                                      double* inBoxCenter){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
158
    //typedef to lighten the writing
159
    typedef FChebSymKernel<double,FChebCell<double,7>, FP2PParticleContainerIndexed<double>, FInterpMatrixKernelR<double>,7> KernelClass;
160
    ChebKernelStruct * newKernel = new ChebKernelStruct();
161
    newKernel->matrix= new FInterpMatrixKernelR<double>();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
162 163 164 165
    int nb_threads = omp_get_max_threads();
    newKernel->kernel = new KernelClass*[nb_threads];
    newKernel->kernel[0]=new KernelClass(inTreeHeight,
                                         inBoxWidth,
166
                                         FPoint<double>(inBoxCenter[0], inBoxCenter[1], inBoxCenter[2]),
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
167 168 169 170 171
                                         newKernel->matrix);

    for(int idThreads=1 ; idThreads<nb_threads ; ++idThreads){
        newKernel->kernel[idThreads] = new KernelClass(*newKernel->kernel[0]);
    }
172 173 174 175 176
    return newKernel;
}

extern "C" void ChebKernelStruct_free(void *inKernel){
    delete reinterpret_cast<ChebKernelStruct *>(inKernel)->matrix;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
177 178 179 180 181
    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;
182 183 184 185
    delete reinterpret_cast<ChebKernelStruct *>(inKernel);
}


186
extern "C" void ChebKernel_P2M(void * leafCell, FSize nbParticles, const FSize *particleIndexes, void * inKernel){
187
    //make temporary array of parts
188
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
189
    tempContainer->reserve(nbParticles);
190
    FPoint<double> pos;
191
    for(int i=0 ; i<nbParticles ; ++i){
192
        pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
193 194
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
195 196 197 198 199 200
        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);
201
    FChebCell<double,7> * realCell = realCellStruct->cell;
202

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
203 204 205
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
210
    realKernel->kernel[id_thread]->P2M(realCell, tempContainer);
211 212 213 214 215 216 217 218
    delete tempContainer;
}

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
219 220
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
221

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
222 223 224
    //Identify thread number
    int id_thread = omp_get_thread_num();

225 226
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
227
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyM2M(childPosition,
228 229
                                                                        childChebCell->getMultipole(0),
                                                                        parentChebCell->getMultipole(0));
230 231 232 233 234 235
}

extern "C" void ChebKernel_M2L(int level, void* targetCell,  void* sourceCell[343], void* inKernel){
    //Get our structures
    ChebCellStruct * targetCellStruct = reinterpret_cast<ChebCellStruct *>(targetCell);
    //get real cheb cell
236
    FChebCell<double,7>* const targetChebCell = targetCellStruct->cell;
237 238

    //copy to an array of FChebCell
239
    const FChebCell<double,7>* arrayOfChebCell[343];
240 241 242 243 244 245 246 247
    for(int i=0; i<343 ; ++i){
        if(sourceCell[i] != nullptr){
            arrayOfChebCell[i] = reinterpret_cast<ChebCellStruct*>(sourceCell[i])->cell;
        }
        else{
            arrayOfChebCell[i] = nullptr;
        }
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
248 249 250 251

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

252 253
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
254
    inKernelStruct->kernel[id_thread]->M2L(targetChebCell,arrayOfChebCell,0,level);
255 256 257
}

extern "C" void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* inKernel){
258
    //Get our structures
259 260 261
    ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
    ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
    //get real cheb cell
262 263
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
264

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
265 266 267
    //Identify thread number
    int id_thread = omp_get_thread_num();

268 269
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
270
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyL2L(childPosition,
271 272
                                                                        parentChebCell->getLocal(0),
                                                                        childChebCell->getLocal(0));
273 274
}

275
extern "C" void ChebKernel_L2P(void* leafCell, FSize nbParticles, const FSize* particleIndexes, void* inKernel){
276
    //Create temporary FSimpleLeaf
277
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
278
    tempContainer->reserve(nbParticles);
279
    FPoint<double> pos;
280
    for(int i=0 ; i<nbParticles ; ++i){
281 282 283 284 285 286
        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);
    }
287 288 289
    //Get our structures
    ChebCellStruct * leafCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    //get real cheb cell
290
    FChebCell<double,7>* leafChebCell = leafCellStruct->cell;
291

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
292 293 294
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
298
    inKernelStruct->kernel[id_thread]->L2P(leafChebCell,tempContainer);
299 300

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

304
    const FVector<FSize>& indexes = tempContainer->getIndexes();
305
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
306 307 308
        forcesToFill[indexes[idxPart]*3+0] += tempContainer->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContainer->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContainer->getForcesZ()[idxPart];
309
        potentialsToFill[indexes[idxPart]] += tempContainer->getPotentials()[idxPart];
310
    }
311 312 313 314 315 316

    delete tempContainer;
    tempContainer=nullptr;
}


317 318
void ChebKernel_P2P(FSize nbParticles, const FSize* particleIndexes,
                    const FSize * sourceParticleIndexes[27],FSize sourceNbPart[27],void* inKernel){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
319

320
    //Create temporary FSimpleLeaf for target
321
    FP2PParticleContainerIndexed<double>* tempContTarget = new FP2PParticleContainerIndexed<double>();
322 323
    tempContTarget->reserve(nbParticles);
    for(int i=0 ; i<nbParticles ; ++i){
324
        FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
325 326
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
327 328 329 330 331
        double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
        tempContTarget->push(pos,particleIndexes[i],Phi);
    }

    //Create 27 FSimpleLeaf for 27 sources
332
    FP2PParticleContainerIndexed<double>* tempContSources[27];
333 334 335
    for(int idSource=0; idSource<27 ; ++idSource){
        if(sourceNbPart[idSource] != 0){
            //Create container
336
            tempContSources[idSource] = new FP2PParticleContainerIndexed<double>();
337 338 339
            //Allocate memory
            tempContSources[idSource]->reserve(sourceNbPart[idSource]);
            //Store a ptr to the indices of that source leaf
340
            const FSize * indSource = sourceParticleIndexes[idSource];
341 342
            //Then, for each part in this source
            for(int i=0 ; i<sourceNbPart[idSource] ; ++i){
343
                FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3  ],
344 345
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+1],
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+2]);
346 347 348 349 350 351 352 353 354 355
                double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[indSource[i]];
                tempContSources[idSource]->push(pos,indSource[i],Phi);
            }
        }
        else{
            tempContSources[idSource] = nullptr;
        }
    }
    //Everything is fine, now, call Chebyshev P2P

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
356 357 358
    //Identify thread number
    int id_thread = omp_get_thread_num();

359 360 361 362 363 364
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
365
    inKernelStruct->kernel[id_thread]->P2P(FTreeCoordinate(coord),tempContTarget,nullptr,tempContSources,0);
366

367
    //get back forces & potentials
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
368
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
369 370
    double * potentialsToFill = reinterpret_cast<UserData *>(inKernel)->potentials[id_thread];

371
    const FVector<FSize>& indexes = tempContTarget->getIndexes();
372
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
373 374 375
        forcesToFill[indexes[idxPart]*3+0] += tempContTarget->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContTarget->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContTarget->getForcesZ()[idxPart];
376
        potentialsToFill[indexes[idxPart]] += tempContTarget->getPotentials()[idxPart];
377 378 379 380 381
    }

    //Note that sources are also modified.
    //get back sources forces
    for(int idSource = 0 ; idSource < 27 ; ++idSource){
382 383
        const FVector<FSize>& indexesSource = tempContSources[idSource]->getIndexes();
        const FSize nbPartInSource = sourceNbPart[idSource];
384
        for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
385 386 387
            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];
388
            potentialsToFill[indexesSource[idxSourcePart]] += tempContSources[idSource]->getPotentials()[idxSourcePart];
389 390 391 392 393 394 395 396 397 398
        }
    }

    //Release memory
    for(int idSource=0; idSource<27 ; ++idSource){
        if(tempContSources[idSource]) delete tempContSources[idSource];
    }
    delete tempContTarget;
}

399
void ChebCell_reset(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell,void * inKernel){
400 401 402 403
    ChebCellStruct *  cellStruct = reinterpret_cast<ChebCellStruct *>(userCell);
    FChebCell<double,7>* chebCell = cellStruct->cell;
    chebCell->resetToInitialState();
}
404 405

#endif