diff --git a/Addons/CKernelApi/Src/CScalfmmApi.h b/Addons/CKernelApi/Src/CScalfmmApi.h index e59af830c8cd58f6210b34b2fe012bb2558ae037..b2e837cb80daf40543212394dd87b591b6aa26e9 100644 --- a/Addons/CKernelApi/Src/CScalfmmApi.h +++ b/Addons/CKernelApi/Src/CScalfmmApi.h @@ -362,7 +362,7 @@ typedef void (*Callback_L2L)(int level, void* parentCell, int childPosition, voi * @param particleIndexes indexes of particles currently computed * @param userData datas specific to the user's kernel */ -typedef void (*Callback_L2P)(void* leafCell, int nbParticles, int* particleIndexes, void* userData); +typedef void (*Callback_L2P)(void* leafCell, int nbParticles,const int* particleIndexes, void* userData); /** * @brief Function to be filled by user's P2P @@ -380,7 +380,7 @@ typedef void (*Callback_P2P)(int nbParticles, const int* particleIndexes, int nb * @param particleIndexes indexes of particles currently computed * @param userData datas specific to the user's kernel */ -typedef void (*Callback_P2PInner)(int nbParticles, int* particleIndexes, void* userData); +typedef void (*Callback_P2PInner)(int nbParticles, const int* particleIndexes, void* userData); @@ -434,6 +434,10 @@ typedef void* (*Callback_init_cell)(int level, long long morton_index, int* tree typedef void (*Callback_free_cell)(void*); +void scalfmm_init_cell(scalfmm_handle Handle, Callback_init_cell user_cell_initializer); + +void scalfmm_free_cell(scalfmm_handle Handle, Callback_free_cell user_cell_deallocator); + ///////////////// Common kernel part : ///////////////// diff --git a/Addons/CKernelApi/Src/FInterEngine.hpp b/Addons/CKernelApi/Src/FInterEngine.hpp index e9a1b66cbfa2e0d298495c3db16b62d161ef8ece..b68c9f1f9832592f41756bbdb4fb518d632f2c34 100644 --- a/Addons/CKernelApi/Src/FInterEngine.hpp +++ b/Addons/CKernelApi/Src/FInterEngine.hpp @@ -67,11 +67,12 @@ public: * @param BoxCenter double[3] coordinate of the center of the * simulation box */ - FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter) : kernel(nullptr), matrix(nullptr), octree(nullptr),arranger(nullptr){ + FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter,scalfmm_kernel_type KernelType) : + kernel(nullptr), matrix(nullptr), octree(nullptr),arranger(nullptr){ octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint(BoxCenter)); this->matrix = new MatrixKernelClass(); this->kernel = new InterKernel(TreeHeight,BoxWidth,FPoint(BoxCenter),matrix); - kernelType = chebyshev; + kernelType = KernelType; arranger = nullptr; } diff --git a/Addons/CKernelApi/Src/FScalFMMEngine.hpp b/Addons/CKernelApi/Src/FScalFMMEngine.hpp index 0c2a51e498690ecec10978d97583455139c9301c..e76ba15f0b9f6e87205703a6562f224670e4152b 100644 --- a/Addons/CKernelApi/Src/FScalFMMEngine.hpp +++ b/Addons/CKernelApi/Src/FScalFMMEngine.hpp @@ -56,7 +56,6 @@ protected: public: FScalFMMEngine() : Algorithm(multi_thread), progress(), nbPart(0){ - //Default behavior in case of out of the box particles is exiting } //First function displayed there are common function for every @@ -197,6 +196,14 @@ public: FAssertLF(0,"No user kernel defined, exiting ...\n"); } + virtual void init_cell(Callback_init_cell user_cell_initializer){ + FAssertLF(0,"No user kernel defined, exiting ...\n"); + } + + virtual void free_cell(Callback_free_cell user_cell_initializer){ + FAssertLF(0,"No user kernel defined, exiting ...\n"); + } + virtual void execute_fmm(){ FAssertLF(0,"No kernel set, cannot execute anything, exiting ...\n"); } @@ -371,6 +378,13 @@ extern "C" void scalfmm_execute_fmm(scalfmm_handle Handle){ ((ScalFmmCoreHandle * ) Handle)->engine->execute_fmm(); } +extern "C" void scalfmm_init_cell(scalfmm_handle Handle, Callback_init_cell user_cell_initializer){ + ((ScalFmmCoreHandle * ) Handle)->engine->init_cell(user_cell_initializer); +} + +extern "C" void scalfmm_free_cell(scalfmm_handle Handle, Callback_free_cell user_cell_deallocator){ + ((ScalFmmCoreHandle * ) Handle)->engine->free_cell(user_cell_deallocator); +} #endif diff --git a/Addons/CKernelApi/Src/FScalfmmApiInit.cpp b/Addons/CKernelApi/Src/FScalfmmApiInit.cpp index 5b7570a0bfae8d6499285f42697f4204ae3dcfe8..03c1154040b87135c898beb8c25fd121cf8415b2 100644 --- a/Addons/CKernelApi/Src/FScalfmmApiInit.cpp +++ b/Addons/CKernelApi/Src/FScalfmmApiInit.cpp @@ -8,6 +8,7 @@ extern "C" { } #include "FInterEngine.hpp" +#include "FUserKernelEngine.hpp" extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* BoxCenter, scalfmm_kernel_type KernelType){ ScalFmmCoreHandle * handle = new ScalFmmCoreHandle(); @@ -26,7 +27,7 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo typedef FInterpMatrixKernelR MatrixKernelClass; typedef FChebSymKernel ChebKernel; - handle->engine = new FInterEngine(TreeHeight,BoxWidth,BoxCenter); + handle->engine = new FInterEngine(TreeHeight,BoxWidth,BoxCenter, KernelType); break; case 2: //TODO typedefs @@ -36,7 +37,7 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo typedef FInterpMatrixKernelR MatrixKernelClass; typedef FUnifKernel UnifKernel; - handle->engine = new FInterEngine(TreeHeight,BoxWidth,BoxCenter); + handle->engine = new FInterEngine(TreeHeight,BoxWidth,BoxCenter,KernelType); break; default: diff --git a/Addons/CKernelApi/Src/FUserKernelEngine.hpp b/Addons/CKernelApi/Src/FUserKernelEngine.hpp new file mode 100644 index 0000000000000000000000000000000000000000..bafd402563b1c2604040c1b4d642baf2c9097b19 --- /dev/null +++ b/Addons/CKernelApi/Src/FUserKernelEngine.hpp @@ -0,0 +1,262 @@ +// =================================================================================== +// Copyright ScalFmm 2014 I +// This software is a computer program whose purpose is to compute the FMM. +// +// This software is governed by the CeCILL-C and LGPL licenses and +// abiding by the rules of distribution of free software. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public and CeCILL-C Licenses for more details. +// "http://www.cecill.info". +// "http://www.gnu.org/licenses". +// =================================================================================== + + +/** + * @file This file contains a class that inherits from FScalFMMEngine, + * and will implement the API functions for a user defined kernel. + */ +#ifndef FUSERKERNELENGINE_HPP +#define FUSERKERNELENGINE_HPP + +#include "FScalFMMEngine.hpp" + + +/** + * @brief CoreCell : Cell used to store User datas + */ +class CoreCell : public FBasicCell { + // Mutable in order to work with the API + mutable void* userData; + +public: + CoreCell() : userData(nullptr) { + } + + ~CoreCell(){ + } + /** + * @brief setContainer store the ptr to the user data inside our + * struct + */ + void setContainer(void* inContainer) const { + userData = inContainer; + } + + /** + * @brief getContainer : return the user datas (in order to give + * it back to the user defined kernel function) + */ + void* getContainer() const { + return userData; + } +}; + + +/** + * This class simply call the function pointers from Scalfmm_Kernel_Descriptor. + * If not pointer is set the calls are skipped. + * The userData is given at any calls. + */ +template< class CellClass, class ContainerClass> +class CoreKernel : public FAbstractKernels { + Scalfmm_Kernel_Descriptor kernel; + void* userData; + +public: + CoreKernel(Scalfmm_Kernel_Descriptor inKernel, void* inUserData) : kernel(inKernel) , userData(inUserData){ + } + + /** Default destructor */ + virtual ~CoreKernel(){ + } + + /** Do nothing */ + virtual void P2M(CellClass* const cell, const ContainerClass* const container) { + if(kernel.p2m) kernel.p2m(cell->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData); + } + + /** Do nothing */ + virtual void M2M(CellClass* const FRestrict cell, const CellClass*const FRestrict *const FRestrict children, const int level) { + if(kernel.m2m){ + for(int idx = 0 ; idx < 8 ; ++idx){ + if( children[idx] ){ + kernel.m2m(level, cell->getContainer(), idx, children[idx]->getContainer(), userData); + } + } + } + } + + /** Do nothing */ + virtual void M2L(CellClass* const FRestrict cell, const CellClass* interactions[], const int , const int level) { + if(kernel.m2l){ + for(int idx = 0 ; idx < 343 ; ++idx){ + if( interactions[idx] ){ + kernel.m2l(level, cell->getContainer(), idx, interactions[idx]->getContainer(), userData); + } + } + } + } + + /** Do nothing */ + virtual void L2L(const CellClass* const FRestrict cell, CellClass* FRestrict *const FRestrict children, const int level) { + if(kernel.l2l){ + for(int idx = 0 ; idx < 8 ; ++idx){ + if( children[idx] ){ + kernel.l2l(level, cell->getContainer(), idx, children[idx]->getContainer(), userData); + } + } + } + } + + /** Do nothing */ + virtual void L2P(const CellClass* const cell, ContainerClass* const container){ + if(kernel.l2p) kernel.l2p(cell->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData); + } + + + /** Do nothing */ + virtual void P2P(const FTreeCoordinate& , + ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/, + ContainerClass* const neighbors[27], const int ){ + if(kernel.p2pinner) kernel.p2pinner(targets->getNbParticles(), targets->getIndexes().data(), userData); + + if(kernel.p2p){ + for(int idx = 0 ; idx < 27 ; ++idx){ + if( neighbors[idx] ){ + kernel.p2p(targets->getNbParticles(), targets->getIndexes().data(), + neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData); + } + } + } + } + + /** Do nothing */ + virtual void P2PRemote(const FTreeCoordinate& , + ContainerClass* const FRestrict , const ContainerClass* const FRestrict , + ContainerClass* const [27], const int ){ + } + +}; + + + + +class FUserKernelEngine : public FScalFMMEngine{ + +private: + + //Typedefs : + typedef FP2PParticleContainerIndexed<> ContainerClass; + typedef FSimpleLeaf LeafClass; + typedef FOctree OctreeClass; + typedef CoreKernel CoreKernelClass; + typedef FP2PLeafInterface LeafInterface; + //For arranger classes + typedef FOctreeArranger ArrangerClass; + typedef FArrangerPeriodic ArrangerClassPeriodic; + + + //Attributes + OctreeClass * octree; + CoreKernelClass * kernel; + ArrangerClass * arranger; + +public: + FUserKernelEngine(int TreeHeight, double BoxWidth , double * BoxCenter, scalfmm_kernel_type KernelType) : + octree(nullptr), kernel(nullptr), arranger(nullptr){ + octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint(BoxCenter)); + kernelType = KernelType; + //Kernel is not set now because the user must provide a + //Scalfmm_Kernel_descriptor + //kernel = new CoreKernelClass(); + } + + + void user_kernel_config( Scalfmm_Kernel_Descriptor userKernel, void * userDatas){ + kernel = new CoreKernelClass(userKernel,userDatas); + } + + + void tree_insert_particles( int NbPositions, double * arrayX, double * arrayY, double * arrayZ){ + for(int idPart = 0; idPartinsert(FPoint(arrayX[idPart],arrayY[idPart],arrayZ[idPart]),idPart); + } + nbPart += NbPositions; + } + + void tree_insert_particles_xyz( int NbPositions, double * XYZ){ + for(int idPart = 0; idPartinsert(FPoint(&XYZ[3*idPart]),idPart); + } + nbPart += NbPositions; + } + + /* + * Call the user allocator on userDatas member field of each cell + */ + void init_cell(Callback_init_cell user_cell_initializer){ + if(user_cell_initializer){ + double boxwidth = octree->getBoxWidth(); + FPoint BoxCenter = octree->getBoxCenter(); + double boxCorner[3]; + boxCorner[0] = BoxCenter.getX() - boxwidth/2.0; + boxCorner[1] = BoxCenter.getY() - boxwidth/2.0; + boxCorner[2] = BoxCenter.getZ() - boxwidth/2.0; + //apply user function on each cell + octree->forEachCellWithLevel([&](CoreCell * currCell,const int currLevel){ + FTreeCoordinate currCoord = currCell->getCoordinate(); + int arrayCoord[3] = {currCoord.getX(),currCoord.getY(),currCoord.getZ()}; + MortonIndex currMorton = currCoord.getMortonIndex(currLevel); + double position[3]; + position[0] = boxCorner[0] + currCoord.getX()*boxwidth/double(1<setContainer(user_cell_initializer(currLevel,currMorton,arrayCoord,position)); + }); + } + } + + void free_cell(Callback_free_cell user_cell_deallocator){ + octree->forEachCell([&](CoreCell * currCell){ + user_cell_deallocator(currCell->getContainer()); + }); + } + + void execute_fmm(){ + //FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n"); + switch(Algorithm){ + case 0: + { + typedef FFmmAlgorithm AlgoClassSeq; + AlgoClassSeq algoSeq(octree,kernel); + algoSeq.execute(); + break; + } + case 1: + { + typedef FFmmAlgorithmThread AlgoClassThread; + AlgoClassThread algoThread(octree,kernel); + algoThread.execute(); + break; + } + case 2: + { + typedef FFmmAlgorithmPeriodic AlgoClassPeriodic; + AlgoClassPeriodic algoPeriod(octree,2); + algoPeriod.setKernel(kernel); + algoPeriod.execute(); + break; + } + default : + std::cout<< "No algorithm found (probably for strange reasons) : "<< Algorithm <<" exiting" << std::endl; + } + + } + +}; + + +#endif //FUSERKERNELENGINE_HPP diff --git a/Addons/CKernelApi/Tests/testUserDefinedKernelApi.c b/Addons/CKernelApi/Tests/testUserDefinedKernelApi.c new file mode 100644 index 0000000000000000000000000000000000000000..ba4e6398d87da45ffa841771ab400a9f4606a2bf --- /dev/null +++ b/Addons/CKernelApi/Tests/testUserDefinedKernelApi.c @@ -0,0 +1,275 @@ +// =================================================================================== +// Logiciel initial: ScalFmm Version 0.5 +// Co-auteurs : Olivier Coulaud, Bérenger Bramas. +// Propriétaires : INRIA. +// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire. +// Initial software: ScalFmm Version 0.5 +// Co-authors: Olivier Coulaud, Bérenger Bramas. +// Owners: INRIA. +// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license. +// =================================================================================== + +#include +#include +#include + +/** + * In this file we implement a example kernel which is simply printing information + * about the cells and particles. + * It is recommanded to compile and execute this code in order to understand the behavior + * of the application. + * We mark some part with "JUST-PUT-HERE:" to give instruction to user to create its own kernel. + **/ + +// Include the FMM API (should be in a different folder for you) +#include "../Src/CScalfmmApi.h" + +// Uncomment the next line to avoid the print in the kernel and verbosity +// #define NOT_TOO_MUCH_VERBOSE +#ifdef NOT_TOO_MUCH_VERBOSE +#define VerbosePrint(X...) +#else +#define VerbosePrint(X...) printf(X) +#endif + +/////////////////////////////////////////////////////////////////////////// +/// Cell struct and functions +/////////////////////////////////////////////////////////////////////////// + +// This represent a cell +struct MyCellDescriptor{ + int level; + long long mortonIndex; + int coord[3]; + double position[3]; + // JUST-PUT-HERE: + // You local and multipole arrays +}; + +// This is our function that init a cell (struct MyCellDescriptor) +void* my_Callback_init_cell(int level, long long mortonIndex, int* coord, double* position){ + VerbosePrint("\tAllocating cell for level %d, morton index %lld, coord %d/%d/%d\n", level, mortonIndex, coord[0], coord[1], coord[2]); + struct MyCellDescriptor* cellData = (struct MyCellDescriptor*)malloc(sizeof(struct MyCellDescriptor)); + memset(cellData, 0, sizeof(struct MyCellDescriptor)); + cellData->level = level; + cellData->mortonIndex = mortonIndex; + memcpy(cellData->coord, coord, sizeof(int)*3); + memcpy(cellData->position, position, sizeof(double)*3); + // JUST-PUT-HERE: + // Fill your structure + return cellData; +} + +// To dealloc a cell +void my_Callback_free_cell(void* cellData){ + free(cellData); +} + + +/////////////////////////////////////////////////////////////////////////// +/// Kernel struct and functions +/////////////////////////////////////////////////////////////////////////// + +// This is the data passed to our kernel +struct MyData { + // We simply strore the number of call the and particle indexes + double* insertedPositions; + int countP2M; + int countM2M; + int countM2L; + int countL2L; + int countL2P; + int countP2PInner; + int countP2P; + // JUST-PUT-HERE: + // everything your kernel will need for its computation + // pre-computed matrices, etc.... +}; + + +// Our P2M +void my_Callback_P2M(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countP2M += 1; + + struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData; + VerbosePrint("Cell morton %lld is doing P2M with %d particles :\n", my_cell->mortonIndex, nbParticlesInLeaf); + int idxPart; + for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3]; + VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart], + position[0], position[1], position[2]); + // JUST-PUT-HERE: + // Your real P2M computation + } +} + +void my_Callback_M2M(int level, void* cellData, int childPosition, void* childData, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countM2M += 1; + + struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData; + struct MyCellDescriptor* my_child = (struct MyCellDescriptor*) childData; + + int childFullPosition[3]; + Scalfmm_utils_parentChildPosition(childPosition, childFullPosition); + + VerbosePrint("Doing a M2M at level %d, between cells %lld and %lld (position %d/%d/%d)\n", + level, my_cell->mortonIndex, my_child->mortonIndex, + childFullPosition[0], childFullPosition[1], childFullPosition[2]); + // JUST-PUT-HERE: your real M2M computation +} + +void my_Callback_M2L(int level, void* cellData, int srcPosition, void* srcData, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countM2L += 1; + + struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData; + struct MyCellDescriptor* my_src_cell = (struct MyCellDescriptor*) srcData; + + int interactionFullPosition[3]; + Scalfmm_utils_interactionPosition(srcPosition, interactionFullPosition); + + VerbosePrint("Doing a M2L at level %d, between cells %lld and %lld (position %d/%d/%d)\n", + level, my_cell->mortonIndex, my_src_cell->mortonIndex, + interactionFullPosition[0], interactionFullPosition[1], interactionFullPosition[2]); + // JUST-PUT-HERE: Your M2L +} + +void my_Callback_L2L(int level, void* cellData, int childPosition, void* childData, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countL2L += 1; + + struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData; + struct MyCellDescriptor* my_child = (struct MyCellDescriptor*) childData; + + int childFullPosition[3]; + Scalfmm_utils_parentChildPosition(childPosition, childFullPosition); + + VerbosePrint("Doing a L2L at level %d, between cells %lld and %lld (position %d/%d/%d)\n", + level, my_cell->mortonIndex, my_child->mortonIndex, + childFullPosition[0], childFullPosition[1], childFullPosition[2]); + // JUST-PUT-HERE: Your L2L +} + +void my_Callback_L2P(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countL2P += 1; + + struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData; + VerbosePrint("Cell morton %lld is doing L2P with %d particles :\n", my_cell->mortonIndex, nbParticlesInLeaf); + int idxPart; + for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3]; + VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart], + position[0], position[1], position[2]); + // JUST-PUT-HERE: Your L2P + } +} + +void my_Callback_P2P(int nbParticlesInLeaf, const int* particleIndexes, int nbParticlesInSrc, const int* particleIndexesSrc, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countP2P += 1; + + VerbosePrint("Doing P2P between two leaves of %d particles and %d particles :\n", nbParticlesInLeaf, nbParticlesInSrc); + int idxPart; + for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3]; + VerbosePrint("\t Target >> particle idx %d position %e/%e%e\n", particleIndexes[idxPart], + position[0], position[1], position[2]); + } + for(idxPart = 0 ; idxPart < nbParticlesInSrc ; ++idxPart){ + double* position = &my_data->insertedPositions[particleIndexesSrc[idxPart]*3]; + VerbosePrint("\t Target >> Src idx %d position %e/%e%e\n", particleIndexesSrc[idxPart], + position[0], position[1], position[2]); + } + + // JUST-PUT-HERE: + // Put one loop into the other to make all particles from source + // interacting with the target particles +} + +void my_Callback_P2PInner(int nbParticlesInLeaf, const int* particleIndexes, void* userData){ + struct MyData* my_data = (struct MyData*)userData; + my_data->countP2PInner += 1; + + VerbosePrint("Doing P2P inside a leaf of %d particles :\n", nbParticlesInLeaf); + int idxPart; + for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3]; + VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart], + position[0], position[1], position[2]); + // JUST-PUT-HERE: another loop to have all particles interacting with + // the other from the same leaf + } +} + +/////////////////////////////////////////////////////////////////////////// +/// Main +/////////////////////////////////////////////////////////////////////////// + +// Simply create particles and try the kernels +int main(int argc, char ** argv){ + // The properties of our tree + int treeHeight = 4; + double boxWidth = 1.0; + double boxCenter[3]; + boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0; + + // Create random particles + int nbParticles = 10; + int particleIndexes[nbParticles]; + double particleXYZ[nbParticles*3]; + { + printf("Creating Particles:\n"); + int idxPart; + for(idxPart = 0 ; idxPart < nbParticles ; ++idxPart){ + particleIndexes[idxPart] = idxPart; + particleXYZ[idxPart*3] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[0]; + particleXYZ[idxPart*3+1] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[1]; + particleXYZ[idxPart*3+2] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[2]; + VerbosePrint("\t %d] %e/%e/%e\n", idxPart, particleXYZ[idxPart*3], particleXYZ[idxPart*3+1], particleXYZ[idxPart*3+2]); + } + } + + // Init the handle + scalfmm_handle handle = scalfmm_init(treeHeight, boxWidth, boxCenter,user_defined_kernel); + // Insert particles + printf("Inserting particles...\n"); + Scalfmm_insert_array_of_particles(handle, nbParticles, particleIndexes, particleXYZ); + // Init cell + printf("Initizalizing the cells:\n"); + Scalfmm_init_cell(handle, my_Callback_init_cell); + + // Init our callback struct + struct User_Scalfmm_Kernel_Descriptor kernel; + kernel.p2m = my_Callback_P2M; + kernel.m2m = my_Callback_M2M; + kernel.m2l = my_Callback_M2L; + kernel.l2l = my_Callback_L2L; + kernel.l2p = my_Callback_L2P; + kernel.p2pinner = my_Callback_P2PInner; + kernel.p2p = my_Callback_P2P; + + // Init the data to pass to all our callbacks + struct MyData my_data; + memset(&my_data, 0, sizeof(struct MyData)); + my_data.insertedPositions = particleXYZ; + + // Execute the FMM + Scalfmm_execute_kernel(handle, kernel, &my_data); + + // Print the results store in our callback + printf("There was %d P2M\n", my_data.countP2M); + printf("There was %d M2M\n", my_data.countM2M); + printf("There was %d M2L\n", my_data.countM2L); + printf("There was %d L2L\n", my_data.countL2L); + printf("There was %d L2P\n", my_data.countL2P); + printf("There was %d P2PInner\n", my_data.countP2PInner); + printf("There was %d P2P\n", my_data.countP2P); + + // Dealloc the handle + Scalfmm_dealloc_handle(handle,my_Callback_free_cell); + + return 0; +}