FScalfmmApiInit.cpp 14.4 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
extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double* BoxCenter, */scalfmm_kernel_type KernelType){
10
    ScalFmmCoreHandle * handle = new ScalFmmCoreHandle();
11
    typedef double FReal;
12 13
    switch(KernelType){
    case 0:
14
        handle->engine = new FUserKernelEngine<FReal>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
15 16 17 18
        break;

    case 1:
        //TODO typedefs
19
        typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
20
        typedef FChebCell<FReal,7>                                                         ChebCell;
21

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

25
        handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
26
        break;
27 28
    // case 2:
    //     //TODO typedefs
29
    //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
30
    //     typedef FUnifCell<7>                                                         UnifCell;
31

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

35 36
    //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
    //     break;
37 38 39 40 41 42 43 44 45

    default:
        std::cout<< "Kernel type unsupported" << std::endl;
        exit(0);
        break;
    }
    return handle;
}

46 47 48 49 50
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
    ((ScalFmmCoreHandle *) handle)->engine->intern_dealloc_handle(userDeallocator);
    delete ((ScalFmmCoreHandle *) handle)->engine ;
    delete (ScalFmmCoreHandle *) handle;
}
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68

/**
 * 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
69
    FChebCell<double,7> * cell;
70 71 72 73 74 75
}ChebCellStruct;


//How to create/destroy cells
extern "C" ChebCellStruct * ChebCellStruct_create(long long int inIndex,int * position){
    ChebCellStruct * newCell = new ChebCellStruct();
76
    newCell->cell = new FChebCell<double,7>();
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
    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 !!!
93 94
    FChebSymKernel<double,FChebCell<double,7>,FP2PParticleContainerIndexed<double>,FInterpMatrixKernelR<double>,7> ** kernel;
    FInterpMatrixKernelR<double> * matrix;
95 96 97 98 99 100
} ChebKernelStruct;

//Kernel functions
extern "C" ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
                                                      double inBoxWidth,
                                                      double* inBoxCenter){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
101
    //typedef to lighten the writing
102
    typedef FChebSymKernel<double,FChebCell<double,7>, FP2PParticleContainerIndexed<double>, FInterpMatrixKernelR<double>,7> KernelClass;
103
    ChebKernelStruct * newKernel = new ChebKernelStruct();
104
    newKernel->matrix= new FInterpMatrixKernelR<double>();
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
105 106 107 108
    int nb_threads = omp_get_max_threads();
    newKernel->kernel = new KernelClass*[nb_threads];
    newKernel->kernel[0]=new KernelClass(inTreeHeight,
                                         inBoxWidth,
109
                                         FPoint<double>(inBoxCenter[0], inBoxCenter[1], inBoxCenter[2]),
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
110 111 112 113 114
                                         newKernel->matrix);

    for(int idThreads=1 ; idThreads<nb_threads ; ++idThreads){
        newKernel->kernel[idThreads] = new KernelClass(*newKernel->kernel[0]);
    }
115 116 117 118 119
    return newKernel;
}

extern "C" void ChebKernelStruct_free(void *inKernel){
    delete reinterpret_cast<ChebKernelStruct *>(inKernel)->matrix;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
120 121 122 123 124
    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;
125 126 127 128
    delete reinterpret_cast<ChebKernelStruct *>(inKernel);
}


129
extern "C" void ChebKernel_P2M(void * leafCell, FSize nbParticles, const int *particleIndexes, void * inKernel){
130
    //make temporary array of parts
131
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
132
    tempContainer->reserve(nbParticles);
133
    FPoint<double> pos;
134
    for(int i=0 ; i<nbParticles ; ++i){
135
        pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
136 137 138 139 140 141 142 143
                     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);
    }
    //get the real cell struct
    ChebCellStruct * realCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
144
    FChebCell<double,7> * realCell = realCellStruct->cell;
145

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
146 147 148
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
153
    realKernel->kernel[id_thread]->P2M(realCell, tempContainer);
154 155 156 157 158 159 160 161
    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
162 163
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
164

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
165 166 167
    //Identify thread number
    int id_thread = omp_get_thread_num();

168 169
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
170
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyM2M(childPosition,
171 172
                                                             childChebCell->getMultipole(0),
                                                             parentChebCell->getMultipole(0));
173 174 175 176 177 178
}

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
179
    FChebCell<double,7>* const targetChebCell = targetCellStruct->cell;
180 181

    //copy to an array of FChebCell
182
    const FChebCell<double,7>* arrayOfChebCell[343];
183 184 185 186 187 188 189 190
    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
191 192 193 194

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

195 196
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
197
    inKernelStruct->kernel[id_thread]->M2L(targetChebCell,arrayOfChebCell,0,level);
198 199 200
}

extern "C" void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* inKernel){
201
   //Get our structures
202 203 204
    ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
    ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
    //get real cheb cell
205 206
    FChebCell<double,7>* parentChebCell = parentCellStruct->cell;
    FChebCell<double,7>* childChebCell = childCellStruct->cell;
207

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
208 209 210
    //Identify thread number
    int id_thread = omp_get_thread_num();

211 212
    //Get the kernel
    ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
213
    inKernelStruct->kernel[id_thread]->getPtrToInterpolator()->applyL2L(childPosition,
214 215
                                                             parentChebCell->getLocal(0),
                                                             childChebCell->getLocal(0));
216 217
}

218
extern "C" void ChebKernel_L2P(void* leafCell, FSize nbParticles, const int* particleIndexes, void* inKernel){
219
    //Create temporary FSimpleLeaf
220
    FP2PParticleContainerIndexed<double>* tempContainer = new FP2PParticleContainerIndexed<double>();
221
    tempContainer->reserve(nbParticles);
222
    FPoint<double> pos;
223
    for(int i=0 ; i<nbParticles ; ++i){
224
    pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
225 226 227 228 229
        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);
}
230 231 232
    //Get our structures
    ChebCellStruct * leafCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
    //get real cheb cell
233
    FChebCell<double,7>* leafChebCell = leafCellStruct->cell;
234

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
235 236 237
    //Identify thread number
    int id_thread = omp_get_thread_num();

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

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
241
    inKernelStruct->kernel[id_thread]->L2P(leafChebCell,tempContainer);
242 243

    //Then retrieve the results
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
244
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
245 246
    const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
247 248 249 250
    forcesToFill[indexes[idxPart]*3+0] += tempContainer->getForcesX()[idxPart];
    forcesToFill[indexes[idxPart]*3+1] += tempContainer->getForcesY()[idxPart];
    forcesToFill[indexes[idxPart]*3+2] += tempContainer->getForcesZ()[idxPart];
}
251 252 253 254 255 256

    delete tempContainer;
    tempContainer=nullptr;
}


257
    void ChebKernel_P2P(FSize nbParticles, const int* particleIndexes,
258
                    const int * sourceParticleIndexes[27],int sourceNbPart[27],void* inKernel){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
259

260
    //Create temporary FSimpleLeaf for target
261
    FP2PParticleContainerIndexed<double>* tempContTarget = new FP2PParticleContainerIndexed<double>();
262 263
    tempContTarget->reserve(nbParticles);
    for(int i=0 ; i<nbParticles ; ++i){
264
        FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3  ],
265 266 267 268 269 270 271
                            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]];
        tempContTarget->push(pos,particleIndexes[i],Phi);
    }

    //Create 27 FSimpleLeaf for 27 sources
272
    FP2PParticleContainerIndexed<double>* tempContSources[27];
273 274 275
    for(int idSource=0; idSource<27 ; ++idSource){
        if(sourceNbPart[idSource] != 0){
            //Create container
276
            tempContSources[idSource] = new FP2PParticleContainerIndexed<double>();
277 278 279 280 281 282
            //Allocate memory
            tempContSources[idSource]->reserve(sourceNbPart[idSource]);
            //Store a ptr to the indices of that source leaf
            const int * indSource = sourceParticleIndexes[idSource];
            //Then, for each part in this source
            for(int i=0 ; i<sourceNbPart[idSource] ; ++i){
283
                FPoint<double> pos = FPoint<double>(reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3  ],
284 285 286 287 288 289 290 291 292 293 294 295
                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+1],
                                    reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+2]);
                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
296 297 298
    //Identify thread number
    int id_thread = omp_get_thread_num();

299 300 301 302 303 304
    //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
305
    inKernelStruct->kernel[id_thread]->P2P(FTreeCoordinate(coord),tempContTarget,nullptr,tempContSources,0);
306 307

    //get back forces
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
308
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
309 310
    const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
311 312 313 314 315 316 317 318
        forcesToFill[indexes[idxPart]*3+0] += tempContTarget->getForcesX()[idxPart];
        forcesToFill[indexes[idxPart]*3+1] += tempContTarget->getForcesY()[idxPart];
        forcesToFill[indexes[idxPart]*3+2] += tempContTarget->getForcesZ()[idxPart];
    }

    //Note that sources are also modified.
    //get back sources forces
    for(int idSource = 0 ; idSource < 27 ; ++idSource){
319
        const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
        const int nbPartInSource = sourceNbPart[idSource];
        for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
            forcesToFill[indexes[idxSourcePart]*3+0] += tempContSources[idSource]->getForcesX()[idxSourcePart];
            forcesToFill[indexes[idxSourcePart]*3+1] += tempContSources[idSource]->getForcesY()[idxSourcePart];
            forcesToFill[indexes[idxSourcePart]*3+2] += tempContSources[idSource]->getForcesZ()[idxSourcePart];
        }
    }

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



#endif