FScalfmmApiInit.cpp 14.6 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 14
    switch(KernelType){
    case 0:
15
        handle->engine = new FUserKernelEngine<FReal>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
16 17 18 19
        break;

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

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

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

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

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

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

47 48 49 50 51
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;
}
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

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


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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    delete tempContainer;
    tempContainer=nullptr;
}


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

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

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

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

    //get back forces
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
309
    double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed[id_thread];
310
    const FVector<FSize>& indexes = tempContTarget->getIndexes();
311
    for(FSize idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
312 313 314 315 316 317 318 319
        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){
320 321
        const FVector<FSize>& indexesSource = tempContSources[idSource]->getIndexes();
        const FSize nbPartInSource = sourceNbPart[idSource];
322
        for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
323 324 325
            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];
326 327 328 329 330 331 332 333 334 335 336 337 338
        }
    }

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



#endif