FScalfmmApiInit.cpp 17.9 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
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>();
12
    typedef double FReal;
13

14
    if(algo == source_target){
15

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

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

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

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

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

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
            //     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;
        }
    }
    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;
64 65
            //typedef FChebCell<FReal,7>                                                   ChebCell;
            typedef FTypedChebCell<FReal,7>                                                   ChebCell;
66 67 68 69 70
            typedef FSimpleLeaf<FReal,ContainerClass>                                         LeafClass;

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

71
            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType,algo);
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
            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;

        }
90
    }
91
   return handle;
92 93
}

94
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
95 96 97
    ((ScalFmmCoreHandle<double> *) handle)->engine->intern_dealloc_handle(userDeallocator);
    delete ((ScalFmmCoreHandle<double> *) handle)->engine ;
    delete (ScalFmmCoreHandle<double> *) handle;
98
}
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116

/**
 * This parts implements all the function defined in FChebInterface.h
 * using the Chebyshev classes
 */
#ifndef CHEBINTERFACE_HPP
#define CHEBINTERFACE_HPP
extern "C" {
#include "Kernels/Chebyshev/FChebInterface.h"
}


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


typedef struct FChebCell_struct{
    //Store what's needed
117
    FChebCell<double,7> * cell;
118 119 120 121 122 123
}ChebCellStruct;


//How to create/destroy cells
extern "C" ChebCellStruct * ChebCellStruct_create(long long int inIndex,int * position){
    ChebCellStruct * newCell = new ChebCellStruct();
124
    newCell->cell = new FChebCell<double,7>();
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
    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 !!!
141 142
    FChebSymKernel<double,FChebCell<double,7>,FP2PParticleContainerIndexed<double>,FInterpMatrixKernelR<double>,7> ** kernel;
    FInterpMatrixKernelR<double> * matrix;
143 144 145 146 147 148
} ChebKernelStruct;

//Kernel functions
extern "C" ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
                                                      double inBoxWidth,
                                                      double* inBoxCenter){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
149
    //typedef to lighten the writing
150
    typedef FChebSymKernel<double,FChebCell<double,7>, FP2PParticleContainerIndexed<double>, FInterpMatrixKernelR<double>,7> KernelClass;
151
    ChebKernelStruct * newKernel = new ChebKernelStruct();
152
    newKernel->matrix= new FInterpMatrixKernelR<double>();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
153 154 155 156
    int nb_threads = omp_get_max_threads();
    newKernel->kernel = new KernelClass*[nb_threads];
    newKernel->kernel[0]=new KernelClass(inTreeHeight,
                                         inBoxWidth,
157
                                         FPoint<double>(inBoxCenter[0], inBoxCenter[1], inBoxCenter[2]),
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
158 159 160 161 162
                                         newKernel->matrix);

    for(int idThreads=1 ; idThreads<nb_threads ; ++idThreads){
        newKernel->kernel[idThreads] = new KernelClass(*newKernel->kernel[0]);
    }
163 164 165 166 167
    return newKernel;
}

extern "C" void ChebKernelStruct_free(void *inKernel){
    delete reinterpret_cast<ChebKernelStruct *>(inKernel)->matrix;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
168 169 170 171 172
    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;
173 174 175 176
    delete reinterpret_cast<ChebKernelStruct *>(inKernel);
}


177
extern "C" void ChebKernel_P2M(void * leafCell, FSize nbParticles, const FSize *particleIndexes, void * inKernel){
178
    //make temporary array of parts
179
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
180
    tempContainer->reserve(nbParticles);
181
    FPoint<double> pos;
182
    for(int i=0 ; i<nbParticles ; ++i){
183
        pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
184 185
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                             reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
186 187 188 189 190 191
        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);
192
    FChebCell<double,7> * realCell = realCellStruct->cell;
193

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
194 195 196
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
201
    realKernel->kernel[id_thread]->P2M(realCell, tempContainer);
202 203 204 205 206 207 208 209
    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
210 211
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
212

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
213 214 215
    //Identify thread number
    int id_thread = omp_get_thread_num();

216 217
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
218
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyM2M(childPosition,
219 220
                                                                        childChebCell->getMultipole(0),
                                                                        parentChebCell->getMultipole(0));
221 222 223 224 225 226
}

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
227
    FChebCell<double,7>* const targetChebCell = targetCellStruct->cell;
228 229

    //copy to an array of FChebCell
230
    const FChebCell<double,7>* arrayOfChebCell[343];
231 232 233 234 235 236 237 238
    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
239 240 241 242

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

243 244
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
245
    inKernelStruct->kernel[id_thread]->M2L(targetChebCell,arrayOfChebCell,0,level);
246 247 248
}

extern "C" void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* inKernel){
249
    //Get our structures
250 251 252
    ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
    ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
    //get real cheb cell
253 254
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
255

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
256 257 258
    //Identify thread number
    int id_thread = omp_get_thread_num();

259 260
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
261
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyL2L(childPosition,
262 263
                                                                        parentChebCell->getLocal(0),
                                                                        childChebCell->getLocal(0));
264 265
}

266
extern "C" void ChebKernel_L2P(void* leafCell, FSize nbParticles, const FSize* particleIndexes, void* inKernel){
267
    //Create temporary FSimpleLeaf
268
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
269
    tempContainer->reserve(nbParticles);
270
    FPoint<double> pos;
271
    for(int i=0 ; i<nbParticles ; ++i){
272 273 274 275 276 277
        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);
    }
278 279 280
    //Get our structures
    ChebCellStruct * leafCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    //get real cheb cell
281
    FChebCell<double,7>* leafChebCell = leafCellStruct->cell;
282

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
283 284 285
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
289
    inKernelStruct->kernel[id_thread]->L2P(leafChebCell,tempContainer);
290 291

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

295
    const FVector<FSize>& indexes = tempContainer->getIndexes();
296
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
297 298 299
        forcesToFill[indexes[idxPart]*3+0] += tempContainer->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContainer->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContainer->getForcesZ()[idxPart];
300
        potentialsToFill[indexes[idxPart]] += tempContainer->getPotentials()[idxPart];
301
    }
302 303 304 305 306 307

    delete tempContainer;
    tempContainer=nullptr;
}


308 309
void ChebKernel_P2P(FSize nbParticles, const FSize* particleIndexes,
                    const FSize * sourceParticleIndexes[27],FSize sourceNbPart[27],void* inKernel){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
310

311
    //Create temporary FSimpleLeaf for target
312
    FP2PParticleContainerIndexed<double>* tempContTarget = new FP2PParticleContainerIndexed<double>();
313 314
    tempContTarget->reserve(nbParticles);
    for(int i=0 ; i<nbParticles ; ++i){
315
        FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
316 317
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
                                            reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
318 319 320 321 322
        double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
        tempContTarget->push(pos,particleIndexes[i],Phi);
    }

    //Create 27 FSimpleLeaf for 27 sources
323
    FP2PParticleContainerIndexed<double>* tempContSources[27];
324 325 326
    for(int idSource=0; idSource<27 ; ++idSource){
        if(sourceNbPart[idSource] != 0){
            //Create container
327
            tempContSources[idSource] = new FP2PParticleContainerIndexed<double>();
328 329 330
            //Allocate memory
            tempContSources[idSource]->reserve(sourceNbPart[idSource]);
            //Store a ptr to the indices of that source leaf
331
            const FSize * indSource = sourceParticleIndexes[idSource];
332 333
            //Then, for each part in this source
            for(int i=0 ; i<sourceNbPart[idSource] ; ++i){
334
                FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3  ],
335 336
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+1],
                                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+2]);
337 338 339 340 341 342 343 344 345 346
                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
347 348 349
    //Identify thread number
    int id_thread = omp_get_thread_num();

350 351 352 353 354 355
    //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
356
    inKernelStruct->kernel[id_thread]->P2P(FTreeCoordinate(coord),tempContTarget,nullptr,tempContSources,0);
357

358
    //get back forces & potentials
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
359
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
360 361
    double * potentialsToFill = reinterpret_cast<UserData *>(inKernel)->potentials[id_thread];

362
    const FVector<FSize>& indexes = tempContTarget->getIndexes();
363
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
364 365 366
        forcesToFill[indexes[idxPart]*3+0] += tempContTarget->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContTarget->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContTarget->getForcesZ()[idxPart];
367
        potentialsToFill[indexes[idxPart]] += tempContTarget->getPotentials()[idxPart];
368 369 370 371 372
    }

    //Note that sources are also modified.
    //get back sources forces
    for(int idSource = 0 ; idSource < 27 ; ++idSource){
373 374
        const FVector<FSize>& indexesSource = tempContSources[idSource]->getIndexes();
        const FSize nbPartInSource = sourceNbPart[idSource];
375
        for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
376 377 378
            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];
379
            potentialsToFill[indexesSource[idxSourcePart]] += tempContSources[idSource]->getPotentials()[idxSourcePart];
380 381 382 383 384 385 386 387 388 389
        }
    }

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

390 391 392 393 394
void ChebCell_reset(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell){
    ChebCellStruct *  cellStruct = reinterpret_cast<ChebCellStruct *>(userCell);
    FChebCell<double,7>* chebCell = cellStruct->cell;
    chebCell->resetToInitialState();
}
395 396

#endif