diff --git a/Addons/FmmApi/Src/FmmApi.h b/Addons/FmmApi/Src/FmmApi.h new file mode 100644 index 0000000000000000000000000000000000000000..200893e18dedd114608d77f49d73917c2f4921d4 --- /dev/null +++ b/Addons/FmmApi/Src/FmmApi.h @@ -0,0 +1,113 @@ +// =================================================================================== +// 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. +// =================================================================================== +#ifndef FMMAPI_H +#define FMMAPI_H + +enum FmmApiErrors { + FMMAPI_NO_ERROR, + FMMAPI_SUPPORTED_PARAMETER, + FMMAPI_UNSUPPORTED_PARAMETER, + FMMAPI_UNKNOWN_PARAMETER +}; + +////////////////////// Opérateurs FMM Core : ////////////////////////// + +enum FmmApiCoreParameters { + FMMCORE_TREE_HEIGHT, // hombre de niveaux de l'arbre (int) + FMMCORE_ROOT_BOX_WIDTH, // taille de la boîte racine (FReal) + FMMCORE_ROOT_BOX_CENTER, // position du centre de la boîte racine (FReal[3]) + FMMCORE_LEAF_BOX_WIDTH, // taille des boîtes feuilles (FReal) + FMMCORE_POINTS_PER_LEAF, // nombre moyen de points par feuille (FReal) + FMMCORE_MPI_COMMUNICATOR, // communicateur MPI (MPI_Comm) + FMMCORE_THREADS_NUMBER, // nombre de threads (int) + FMMCORE_THREAD_ID, // id du thread (int) + FMMCORE_RHS_NUMBER, // nombre de seconds membres (int) + //paramètres en lecture seule : + FMMCORE_HANDLES_P2P // renvoie 0 ou 1 pour dire si le FmmCore gère ou pas le P2P. +}; + +int FmmCore_init(void **fmmCore) ; /*alloue et initialise le FmmCore*/ +int FmmCore_free(void *fmmCore) ; /*libère le FmmCore*/ +int FmmCore_isParameterUsed(void */*fmmCore*/, int *name, int *flag); +int FmmCore_setParameter(void *fmmCore, int *name, void*value); +int FmmCore_setParameter(void *fmmCore, int name, void*value); +int FmmCore_getParameter(void *fmmCore, int *name, void*value); +int FmmCore_getParameter(void *fmmCore, int name, void*value); + + +int FmmCore_getRadius(void*fmmCore, void *boxId, FReal *radius);/*Renvoie le rayon de la boîte*/ +int FmmCore_getCentre(void*/*fmmCore*/, void *boxId, FReal **centre); /*Renvoie dans le FReal[3] centre les coordonnées du centre*/ +int FmmCore_getLevel(void*/*fmmCore*/, void *boxId, int *level); /*Renvoie dans level le niveau de la boîte boxId dans son arbre*/ +int FmmCore_getMultipoleArray(void* /*fmmCore*/, void *boxId, void **F); /*Renvoie dans F l'adresse où stocker l'expansion multipôle associée à la boîte boxId*/ +int FmmCore_getLocalArray(void* /*fmmCore*/, void *boxId, void **F); /*Renvoie dans F l'adresse où stocker l'expansion locale associée à la boîte boxId*/ +int FmmCore_getCoord(void*/*fmmCore*/, void *boxId, int *coord); /*Renvoie dans coord la position dans l'arbre*/ + + +/* Données potentiel/champ */ +int FmmCore_getSource(void* /*fmmCore*/, void *boxId, FReal** position, void** potential, int *number); /* Appelé par P2P et P2M pour obtenir le nombre, la position et le potentiel des sources. +Les différents tableaux sont (éventuellement) alloués par le FmmCore. */ +int FmmCore_releaseSource(void*fmmCore, void *boxId, void* potential, FReal* position); /* si le core veut libérer ces tableaux potentiel et position.*/ +int FmmCore_getTargetPoints(void* /*fmmCore*/, void *boxId, FReal** position, int *number) ; /* : Appelé par P2P et L2P pour obtenir le nombre et la position des points cibles.*/ +int FmmCore_releaseTargetPoints(void*fmmCore, void *boxId, FReal* position); /* si le core veut libérer ce tableau "position".*/ +int FmmCore_getTargetField(void* /*fmmCore*/, void *boxId, FReal* F); /* obtient dans un tableau F alloué/libéré par L2P/P2P les valeurs des champs aux points cible +(pour le cas où P2P et L2P doivent sommer leurs résultats).*/ +int FmmCore_setTargetField(void* /*fmmCore*/, void *boxId, FReal* F); /* transmets au FmmCore dans F les valeurs des champs aux points cibles mis à jour.*/ + +/* Entrée/sortie principale */ +int FmmCore_setKernelData(void *fmmCore, void *fmmKernel); /* stocke l'identifiant du FmmKernel dans le FmmCore. Cela permet par la suite aux +opérateurs FMM d'accéder au FmmKernel, et donc d'accéder aux données spécifiques au kernel +(p.ex. fréquence dans le cas Helmholtz, …)*/ +int FmmCore_getKernelData(void*fmmCore, void **fmmKernel); /* récupère l'identifiant du FmmKernel. */ +int FmmCore_setPositions(void *fmmCore, int *nb, FReal *position) ; /* transmet au FmmCore les potentiels associés aux points sources. +Le tableau potential est alloué et libéré par la routine appelant, le FmmCore doit donc en faire une copie.*/ +int FmmCore_setPotentials(void *fmmCore, void *potentials); /* transmet au FmmCore les potentiels associés aux points sources. +Le tableau potential est alloué et libéré par la routine appelant, le FmmCore doit donc en faire une copie.*/ +int FmmCore_doComputation(void *fmmCore) ; /* réalise le produit multipôle. */ +/* !!! Warning use *filed and not **field */ +int FmmCore_getField(void *fmmCore, void *fields) ;/* récupère après le produit multipôle la valeur des champs en chaque point. +Le tableau field doit être alloué et libéré par la routine appelante. */ + +////////////////////// Opérateurs FMM Kernel : ////////////////////////// + +enum FmmApiKernelParameters { + FMMKERNEL_ACCURACY, // précision demandée à la FMM : 1e-3, 1e-6, ... (FReal) + FMMKERNEL_POTENTIAL_DATA_SIZE, // taille en octet de la donnée "potentiel" pour 1 point source (int) + FMMKERNEL_FIELD_DATA_SIZE, // taille en octet de la donnée "field" pour 1 point cible (int) + //paramètres en lecture seule : + FMMKERNEL_HANDLES_P2P // renvoie 0 ou 1 pour dire si le FmmKernel gère ou pas le P2P. +}; + +/******* Allocation : ******/ +int FmmKernel_init(void *fmmCore, void **fmmKernel);/* : alloue et initialise le FmmKernel */ +int FmmKernel_free(void *fmmKernel); /* libére le FmmKernel */ + +/******* Configuration : ***/ + +int FmmKernel_isParameterUsed(void * /*fmm*/, int *name, int *flag); +int FmmKernel_setParameter(void *fmmKernel, int *name, void*value); +int FmmKernel_setParameter(void *fmmKernel, int name, void*value); +int FmmKernel_getParameter(void *fmmKernel, int *name, void*value); +int FmmKernel_getParameter(void *fmmKernel, int name, void*value); + +/****** Données FMM : *****/ +int FmmKernel_getMultipoleArraySize(void *fmmCore, int *size); /* Renvoie dans size la taille (en octets) de l'expansion multipôle associée à la boîte boxId */ +int FmmKernel_getLocalArraySize(void *fmmCore, int *size); /* Renvoie dans size la taille (en octets) de l'expansion locale associée à la boîte boxId*/ + +/******* Opérateurs FMM : **/ +int FmmKernel_P2M(void *fmmCore, void* boxId); +int FmmKernel_L2P(void *fmmCore, void* boxId); +int FmmKernel_M2M(void *fmmCore, void *boxIdFather, void *boxIdSon); +int FmmKernel_L2L(void *fmmCore, void *boxIdFather, void *boxIdSon); +int FmmKernel_M2L(void *fmmCore, void *boxIdSrc, void *boxIdDest); +int FmmKernel_P2P(void *fmmCore, void *boxIdSrc, void *boxIdDest); /* pas mutuel, i.e. on fait seulement dans 1 sens. */ + + +#endif // FMMAPI_H diff --git a/Addons/FmmApi/Src/ScalfmmApiCore.cpp b/Addons/FmmApi/Src/ScalfmmApiCore.cpp new file mode 100644 index 0000000000000000000000000000000000000000..673c67ea378db66d99ac6c1e0dab245cc3673a96 --- /dev/null +++ b/Addons/FmmApi/Src/ScalfmmApiCore.cpp @@ -0,0 +1,689 @@ +// =================================================================================== +// 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 "../../Src/Utils/FParameters.hpp" +#include "../../Src/Utils/FTic.hpp" + +#include "../../Src/Containers/FOctree.hpp" +#include "../../Src/Containers/FVector.hpp" + +#include "../../Src/Components/FSimpleLeaf.hpp" +#include "../../Src/Components/FBasicParticle.hpp" + +#include "../../Src/Utils/FPoint.hpp" + +#include "../../Src/Core/FFmmAlgorithm.hpp" +#include "../../Src/Core/FFmmAlgorithmThread.hpp" +#include "../../Src/Core/FFmmAlgorithmTask.hpp" + +#include "../../Src/Components/FBasicKernels.hpp" + +#ifdef SCALFMM_USE_MPI +#include "../../Src/Utils/FMpi.hpp" +#endif + +#include "FmmApi.h" + + +template <class ParticleClass> +class CoreParticle : public ParticleClass { + int index; +public: + CoreParticle() : index(0) { + } + void setIndex(const int inIndex){ + index = inIndex; + } + int getIndex() const { + return index; + } +}; + +template <class ContainerClass> +class CoreCell : public FBasicCell { + char* multipole; + char* local; + int level; + FReal positions[3]; + mutable ContainerClass* container; + +public: + CoreCell() : multipole(0), local(0), level(0), container(0) { + } + void createArrays(const int multipoleSize, const int localSize){ + multipole = new char[multipoleSize]; + local = new char[localSize]; + + memset(multipole, 0, multipoleSize); + memset(local, 0, localSize); + } + ~CoreCell(){ + delete[] multipole; + delete[] local; + } + + const void* getMultipole() const{ + return multipole; + } + + const void* getLocal() const{ + return local; + } + + void* getMultipole(){ + return multipole; + } + + void* getLocal(){ + return local; + } + + void setLevel(const int inLevel){ + level = inLevel; + } + + int getLevel() const { + return level; + } + + void setPosition(const FReal inPositions[3]){ + positions[0] = inPositions[0]; + positions[1] = inPositions[1]; + positions[2] = inPositions[2]; + } + + const FReal* getPosition() const { + return positions; + } + + void setContainer(ContainerClass* inContainer) const { + container = inContainer; + } + + ContainerClass* getContainer() const { + return container; + } +}; + + + +template< class ParticleClass, class CellClass, class ContainerClass> +class CoreKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> { + void* fmmCore; + +public: + CoreKernel(void* inFmmCore): fmmCore(inFmmCore){ + } + + /** Default destructor */ + virtual ~CoreKernel(){ + } + + /** Do nothing */ + virtual void P2M(CellClass* const cell, const ContainerClass* const container) { + cell->setContainer(const_cast<ContainerClass*>(container)); + FmmKernel_P2M(fmmCore,(void*) cell); + } + + /** Do nothing */ + virtual void M2M(CellClass* const FRestrict cell, const CellClass*const FRestrict *const FRestrict children, const int ) { + for(int idx = 0 ; idx < 8 ; ++idx){ + if( children[idx] ){ + FmmKernel_M2M(fmmCore, (void*)cell, (void*)children[idx]); + } + } + } + + /** Do nothing */ + virtual void M2L(CellClass* const FRestrict cell, const CellClass* interactions[], const int , const int ) { + for(int idx = 0 ; idx < 343 ; ++idx){ + if( interactions[idx] ){ + FmmKernel_M2L(fmmCore, (void*)cell, (void*)interactions[idx]); + } + } + } + + /** Do nothing */ + virtual void L2L(const CellClass* const FRestrict cell, CellClass* FRestrict *const FRestrict children, const int ) { + for(int idx = 0 ; idx < 8 ; ++idx){ + if( children[idx] ){ + FmmKernel_L2L(fmmCore, (void*)cell, (void*)children[idx]); + } + } + } + + /** Do nothing */ + virtual void L2P(const CellClass* const cell, ContainerClass* const container){ + cell->setContainer((ContainerClass*)container); + FmmKernel_L2P(fmmCore, (void*)cell); + } + + + /** Do nothing */ + virtual void P2P(const FTreeCoordinate& , + ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/, + ContainerClass* const neighbors[27], const int ){ + CellClass* cell = (CellClass*)targets->getParentCell(); + FmmKernel_P2P(fmmCore, cell, cell); + + for(int idx = 0 ; idx < 27 ; ++idx){ + if( neighbors[idx] ){ + FmmKernel_P2P(fmmCore, (CellClass*)neighbors[idx]->getParentCell(), cell); + } + } + } + + /** Do nothing */ + virtual void P2PRemote(const FTreeCoordinate& , + ContainerClass* const FRestrict , const ContainerClass* const FRestrict , + ContainerClass* const [27], const int ){ + printf("Error remote not implemented!!!!\n"); + } + +}; + + +template <class ObjectClass> +class CoreVector { + mutable void* parentCell; + void* fields; + int sizeOfField; + void* potentials; + int sizeOfPential; + FVector<FReal> positions; + FVector<int> indexes; +public: + CoreVector() : parentCell(0), fields(0), sizeOfField(0), + potentials(0), sizeOfPential(0){ + } + + ~CoreVector(){ + delete[] (char*)fields; + delete[] (char*)potentials; + } + + void setParentCell(void* inCell) const{ + parentCell = inCell; + } + + void* getParentCell() const{ + return parentCell; + } + + void allocateFields(const int memSizePerParticles){ + if(fields) delete[](char*)fields; + fields = new char[getNbParticles() * memSizePerParticles]; + sizeOfField = memSizePerParticles; + memset(fields, 0, getNbParticles() * memSizePerParticles); + } + + void* getFields() const{ + return fields; + } + + int getSizeOfField() const{ + return sizeOfField; + } + + void allocatePotentials(const int memSizePerParticles){ + if(potentials) delete[](char*)potentials; + potentials = new char[getNbParticles() * memSizePerParticles]; + sizeOfPential = memSizePerParticles; + memset(potentials, 0, getNbParticles() * memSizePerParticles); + } + + void* getPotentials() const{ + return potentials; + } + + int getSizeOfPotential() const{ + return sizeOfPential; + } + + const int* getIndexes() const { + return indexes.data(); + } + + const FReal* getPositions() const { + return positions.data(); + } + + int getNbParticles() { + return indexes.getSize(); + } + + void push(const ObjectClass& particle){ + positions.push(particle.getPosition().getX()); + positions.push(particle.getPosition().getY()); + positions.push(particle.getPosition().getZ()); + indexes.push(particle.getIndex()); + } +}; + +typedef CoreParticle<FBasicParticle> CoreParticleClass; +typedef CoreVector<CoreParticleClass> CoreContainerClass; + +typedef CoreCell<CoreContainerClass> CoreCellClass; +typedef FSimpleLeaf<CoreParticleClass, CoreContainerClass > LeafClass; +typedef FOctree<CoreParticleClass, CoreCellClass, CoreContainerClass , LeafClass > OctreeClass; +typedef CoreKernel<CoreParticleClass, CoreCellClass, CoreContainerClass> CoreKernelClass; + +typedef FFmmAlgorithm<OctreeClass, CoreParticleClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClass; +typedef FFmmAlgorithmThread<OctreeClass, CoreParticleClass, CoreCellClass, CoreContainerClass, CoreKernelClass, LeafClass > FmmClassThread; + + +struct ScalFmmCoreHandle { + struct ScalFmmCoreConfig { + //paramètres en lecture/écriture : + int treeHeight; // hombre de niveaux de l'arbre (int) + FReal boxWidth; // taille de la boîte racine (FReal) + FReal boxCenter[3]; // position du centre de la boîte racine (FReal[3]) +#ifdef SCALFMM_USE_MPI + MPI_Comm mpiCom; // communicateur MPI (MPI_Comm) +#endif + int nbThreads; // nombre de threads (int) + int rhsNumber; // nombre de seconds membres (int) + }; + + ScalFmmCoreConfig config; + OctreeClass* octree; + void *kernelHandle; +}; + +int FmmCore_init(void **fmmCore) { + ScalFmmCoreHandle* corehandle = new ScalFmmCoreHandle; + memset(corehandle, 0, sizeof(corehandle)); + + corehandle->config.nbThreads = omp_get_max_threads(); + + *fmmCore = corehandle; + + return FMMAPI_NO_ERROR; +} /*alloue et initialise le FmmCore*/ + +// !!! Added !!! +int FmmCore_finishInit(void *fmmCore) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*) fmmCore; + + corehandle->octree = new OctreeClass(corehandle->config.treeHeight, FMath::Min(3,corehandle->config.treeHeight-1), + corehandle->config.boxWidth, FPoint(corehandle->config.boxCenter)); + + if( corehandle->config.nbThreads != 0){ + omp_set_num_threads(corehandle->config.nbThreads); + } + + return FMMAPI_NO_ERROR; +} /*alloue et initialise le FmmCore*/ + +int FmmCore_free(void *fmmCore) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + if(corehandle->octree) delete corehandle->octree; + delete corehandle; + return FMMAPI_NO_ERROR; +} /*libère le FmmCore*/ + +int FmmCore_isParameterUsed(void */*fmmCore*/, int *name, int *flag){ + switch( *name ){ + case FMMCORE_ROOT_BOX_WIDTH : + case FMMCORE_ROOT_BOX_CENTER : + case FMMCORE_TREE_HEIGHT : +#ifdef SCALFMM_USE_MPI + case FMMCORE_MPI_COMMUNICATOR: +#endif + case FMMCORE_THREADS_NUMBER: + case FMMCORE_THREAD_ID: + case FMMCORE_RHS_NUMBER: + case FMMCORE_HANDLES_P2P: + *flag = FMMAPI_SUPPORTED_PARAMETER; + break; + case FMMCORE_LEAF_BOX_WIDTH : + case FMMCORE_POINTS_PER_LEAF : + *flag = FMMAPI_UNSUPPORTED_PARAMETER; + break; + default: + *flag = FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmCore_setParameter(void *fmmCore, int *name, void*value){ + int flag; + + FmmCore_isParameterUsed(fmmCore, name, &flag); + if( flag != FMMAPI_SUPPORTED_PARAMETER){ + return flag; + } + + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + switch( *name ){ + case FMMCORE_TREE_HEIGHT : + corehandle->config.treeHeight = *(int*)value ; + break; + case FMMCORE_ROOT_BOX_WIDTH : + corehandle->config.boxWidth = *(FReal*)value; + break; + case FMMCORE_ROOT_BOX_CENTER : + memcpy(corehandle->config.boxCenter, value, sizeof(FReal)*3); + break; +#ifdef SCALFMM_USE_MPI + case FMMCORE_MPI_COMMUNICATOR: + corehandle->config.mpiCom = *(MPI_Comm*)value; + break; +#endif + case FMMCORE_THREADS_NUMBER: + corehandle->config.nbThreads = *(int*)value; + break; + case FMMCORE_RHS_NUMBER: + corehandle->config.rhsNumber = *(int*)value; + break; + case FMMCORE_HANDLES_P2P: + case FMMCORE_THREAD_ID: + return FMMAPI_UNSUPPORTED_PARAMETER; + default: + return FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmCore_setParameter(void *fmmCore, int name, void*value){ + return FmmCore_setParameter(fmmCore, &name, value); +} + + +int FmmCore_getParameter(void *fmmCore, int *name, void*value){ + int flag; + + FmmCore_isParameterUsed(fmmCore, name, &flag); + if( flag != FMMAPI_SUPPORTED_PARAMETER){ + return flag; + } + + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + switch( *name ){ + case FMMCORE_TREE_HEIGHT : + *(int*)value = corehandle->config.treeHeight; + break; + case FMMCORE_ROOT_BOX_WIDTH : + *(FReal*)value = corehandle->config.boxWidth; + break; + case FMMCORE_ROOT_BOX_CENTER : + memcpy(value,corehandle->config.boxCenter, sizeof(FReal)*3); + break; +#ifdef SCALFMM_USE_MPI + case FMMCORE_MPI_COMMUNICATOR: + *(MPI_Comm*)value = corehandle->config.mpiCom; + break; +#endif + case FMMCORE_THREADS_NUMBER: + *(int*)value = corehandle->config.nbThreads; + break; + case FMMCORE_THREAD_ID: + *(int*)value = omp_get_thread_num(); + break; + case FMMCORE_RHS_NUMBER: + *(int*)value = corehandle->config.rhsNumber; + break; + case FMMCORE_HANDLES_P2P: + *(int*)value = true; + break; + default: + return FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmCore_getParameter(void *fmmCore, int name, void*value){ + return FmmCore_getParameter(fmmCore, &name, value); +} + + +int FmmCore_getRadius(void*fmmCore, void *boxId, FReal *radius) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + *radius = corehandle->octree->getBoxWidth() / FReal(1<<boxhandle->getLevel()); + return FMMAPI_NO_ERROR; +} /*Renvoie le rayon de la boîte*/ + +int FmmCore_getCentre(void*/*fmmCore*/, void *boxId, FReal **centre) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + + (*centre)[0] = boxhandle->getPosition()[0]; + (*centre)[1] = boxhandle->getPosition()[1]; + (*centre)[2] = boxhandle->getPosition()[2]; + + return FMMAPI_NO_ERROR; +} /*Renvoie dans le FReal[3] centre les coordonnées du centre*/ + +int FmmCore_getLevel(void*/*fmmCore*/, void *boxId, int *level) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + *level = boxhandle->getLevel(); + return FMMAPI_NO_ERROR; +} /*Renvoie dans level le niveau de la boîte boxId dans son arbre*/ + +int FmmCore_getMultipoleArray(void* /*fmmCore*/, void *boxId, void **F) { + CoreCellClass* cell = (CoreCellClass*)boxId; + *F = cell->getMultipole(); + return FMMAPI_NO_ERROR; +} /*Renvoie dans F l'adresse où stocker l'expansion multipôle associée à la boîte boxId*/ + +int FmmCore_getLocalArray(void* /*fmmCore*/, void *boxId, void **F) { + CoreCellClass* cell = (CoreCellClass*)boxId; + *F = cell->getLocal(); + return FMMAPI_NO_ERROR; +} /*Renvoie dans F l'adresse où stocker l'expansion locale associée à la boîte boxId*/ + +int FmmCore_getCoord(void*/*fmmCore*/, void *boxId, int *coord) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + coord[0] = boxhandle->getCoordinate().getX(); + coord[1] = boxhandle->getCoordinate().getY(); + coord[2] = boxhandle->getCoordinate().getZ(); + return FMMAPI_NO_ERROR; +} /*Renvoie dans coord la position dans l'arbre*/ + + +/* Données potentiel/champ */ +int FmmCore_getSource(void* /*fmmCore*/, void *boxId, FReal** position, void** potential, int *number) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + CoreContainerClass* sources = boxhandle->getContainer(); + + *number = sources->getNbParticles(); + *position = const_cast<FReal*>(sources->getPositions()); + *potential = sources->getPotentials(); + + return FMMAPI_NO_ERROR; +} /* Appelé par P2P et P2M pour obtenir le nombre, la position et le potentiel des sources. +Les différents tableaux sont (éventuellement) alloués par le FmmCore. */ + +int FmmCore_releaseSource(void*fmmCore, void *boxId, void* potential, FReal* position) { + return FMMAPI_NO_ERROR; +} /* si le core veut libérer ces tableaux potentiel et position.*/ + +int FmmCore_getTargetPoints(void* /*fmmCore*/, void *boxId, FReal** position, int *number) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + CoreContainerClass* targets = boxhandle->getContainer(); + + *number = targets->getNbParticles(); + *position = const_cast<FReal*>(targets->getPositions()); + + return FMMAPI_NO_ERROR; +} /* : Appelé par P2P et L2P pour obtenir le nombre et la position des points cibles.*/ + +int FmmCore_releaseTargetPoints(void*fmmCore, void *boxId, FReal* position) { + return FMMAPI_NO_ERROR; +} /* si le core veut libérer ce tableau "position".*/ + +int FmmCore_getTargetField(void* /*fmmCore*/, void *boxId, FReal* F) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + CoreContainerClass* targets = boxhandle->getContainer(); + + memcpy(F, targets->getFields(), targets->getNbParticles() * targets->getSizeOfField()); + + return FMMAPI_NO_ERROR; +} /* obtient dans un tableau F alloué/libéré par L2P/P2P les valeurs des champs aux points cible +(pour le cas où P2P et L2P doivent sommer leurs résultats).*/ + +int FmmCore_setTargetField(void* /*fmmCore*/, void *boxId, FReal* F) { + CoreCellClass* boxhandle = (CoreCellClass*)boxId; + CoreContainerClass* targets = boxhandle->getContainer(); + + memcpy(targets->getFields(), F, targets->getNbParticles() * targets->getSizeOfField()); + + return FMMAPI_NO_ERROR; +} /* transmets au FmmCore dans F les valeurs des champs aux points cibles mis à jour.*/ + + + +/* Entrée/sortie principale */ +int FmmCore_setKernelData(void *fmmCore, void *fmmKernel) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + corehandle->kernelHandle = fmmKernel; + return FMMAPI_NO_ERROR; +} /* stocke l'identifiant du FmmKernel dans le FmmCore. Cela permet par la suite aux +opérateurs FMM d'accéder au FmmKernel, et donc d'accéder aux données spécifiques au kernel +(p.ex. fréquence dans le cas Helmholtz, …)*/ + +int FmmCore_getKernelData(void*fmmCore, void **fmmKernel) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + *fmmKernel = corehandle->kernelHandle; + return FMMAPI_NO_ERROR; +} /* récupère l'identifiant du FmmKernel. */ + + +int FmmCore_setPositions(void *fmmCore, int *nb, FReal *position) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + OctreeClass* octree = corehandle->octree; + + CoreParticleClass part; + for(int idxPart = 0 ; idxPart < (*nb) ; ++idxPart){ + FReal* pos = &position[idxPart * 3]; + part.setPosition(pos[0], pos[1], pos[2]); + part.setIndex(idxPart); + octree->insert(part); + } + + return FMMAPI_NO_ERROR; +} /* transmet au FmmCore les potentiels associés aux points sources. +Le tableau potential est alloué et libéré par la routine appelant, le FmmCore doit donc en faire une copie.*/ + + +int FmmCore_setPotentials(void *fmmCore, void *potentials) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + OctreeClass* octree = corehandle->octree; + + void* fmmKernel; + FmmCore_getKernelData(fmmCore, &fmmKernel); + + int sizeOfPotentials; + FmmKernel_getParameter(fmmKernel, FMMKERNEL_POTENTIAL_DATA_SIZE, &sizeOfPotentials); + + typename OctreeClass::Iterator octreeIterator(octree); + octreeIterator.gotoBottomLeft(); + do{ + CoreContainerClass* particles = octreeIterator.getCurrentLeaf()->getSrc(); + particles->allocatePotentials(sizeOfPotentials); + + for(int idx = 0 ; idx < particles->getNbParticles() ; ++idx){ + memcpy(((char*)particles->getPotentials())+idx*sizeOfPotentials, + ((char*)potentials)+particles->getIndexes()[idx]*sizeOfPotentials, + sizeOfPotentials); + } + + } while( octreeIterator.moveRight() ); + + return FMMAPI_NO_ERROR; +} /* transmet au FmmCore les potentiels associés aux points sources. +Le tableau potential est alloué et libéré par la routine appelant, le FmmCore doit donc en faire une copie.*/ + +int FmmCore_doComputation(void *fmmCore) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + + { // Ceck if there is number of NbPart summed at level 1 + FPoint corner = corehandle->octree->getBoxCenter(); + corner -= (corehandle->octree->getBoxWidth()/2); + + FReal boxWidth = corehandle->octree->getBoxWidth()/( 1<< (corehandle->config.treeHeight-1) ); + + typename OctreeClass::Iterator octreeIterator(corehandle->octree); + octreeIterator.gotoBottomLeft(); + for(int idxLevel = corehandle->config.treeHeight - 1 ; idxLevel > 1 ; --idxLevel ){ + int multipoleSize; + int localSize; + + FmmKernel_getMultipoleArraySize(fmmCore, &multipoleSize); + FmmKernel_getLocalArraySize(fmmCore, &localSize); + + do{ + octreeIterator.getCurrentCell()->createArrays(multipoleSize, localSize); + octreeIterator.getCurrentCell()->setLevel(idxLevel); + const FTreeCoordinate coord = octreeIterator.getCurrentCell()->getCoordinate(); + FPoint position( coord.getX()*boxWidth + boxWidth/2.0 + corner.getX(), + coord.getY()*boxWidth + boxWidth/2.0 + corner.getY(), + coord.getZ()*boxWidth + boxWidth/2.0 + corner.getZ()); + octreeIterator.getCurrentCell()->setPosition(position.getDataValue()); + } while( octreeIterator.moveRight() ); + + octreeIterator.moveUp(); + octreeIterator.gotoLeft(); + + boxWidth *= FReal(2.0); + } + } + { // Ceck if there is number of NbPart summed at level 1 + void* fmmKernel; + FmmCore_getKernelData(fmmCore, &fmmKernel); + + int sizeOfFields; + FmmKernel_getParameter(fmmKernel, FMMKERNEL_FIELD_DATA_SIZE, &sizeOfFields); + + typename OctreeClass::Iterator octreeIterator(corehandle->octree); + octreeIterator.gotoBottomLeft(); + do{ + octreeIterator.getCurrentLeaf()->getTargets()->setParentCell(octreeIterator.getCurrentCell()); + octreeIterator.getCurrentLeaf()->getTargets()->allocateFields(sizeOfFields); + } while( octreeIterator.moveRight() ); + } + + if( corehandle->config.nbThreads <= 1){ + CoreKernelClass kernels(fmmCore ); + FmmClass algo(corehandle->octree,&kernels); + algo.execute(); + } + else{ + CoreKernelClass kernels(fmmCore ); + FmmClassThread algo(corehandle->octree,&kernels); + algo.execute(); + } + + return FMMAPI_NO_ERROR; +} /* réalise le produit multipôle. */ + +/* !!! Warning use *filed and not **field */ +int FmmCore_getField(void *fmmCore, void *fields) { + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + + typename OctreeClass::Iterator octreeIterator(corehandle->octree); + octreeIterator.gotoBottomLeft(); + do{ + CoreContainerClass* particles = octreeIterator.getCurrentLeaf()->getTargets(); + for(int idx = 0 ; idx < particles->getNbParticles() ; ++idx){ + memcpy(((char*)fields)+particles->getIndexes()[idx]*particles->getSizeOfField(), + ((char*)particles->getFields())+idx*particles->getSizeOfField(), + particles->getSizeOfField()); + } + } while(octreeIterator.moveRight()); + + return FMMAPI_NO_ERROR; +} /* récupère après le produit multipôle la valeur des champs en chaque point. +Le tableau field doit être alloué et libéré par la routine appelante. */ + + diff --git a/Addons/FmmApi/Src/ScalfmmApiRotation.cpp b/Addons/FmmApi/Src/ScalfmmApiRotation.cpp new file mode 100644 index 0000000000000000000000000000000000000000..101921bcfbe66060ce332130a90263763ac18729 --- /dev/null +++ b/Addons/FmmApi/Src/ScalfmmApiRotation.cpp @@ -0,0 +1,472 @@ +// =================================================================================== +// 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 "../../Src/Containers/FVector.hpp" + +#include "../../Src/Components/FSimpleLeaf.hpp" + +#include "../../Src/Utils/FPoint.hpp" + +#include "../../Src/Kernels/Rotation/FRotationCell.hpp" +#include "../../Src/Kernels/Rotation/FRotationParticle.hpp" +#include "../../Src/Kernels/Rotation/FRotationKernel.hpp" + +#include "FmmApi.h" + +////////////////////// Opérateurs FMM Kernel : ////////////////////////// + +class KernelCell : public FBasicCell { + FComplexe* multipole; + FComplexe* local; +public: + KernelCell() : multipole(0), local(0){ + } + void attachArrays(FComplexe inMultipole[], FComplexe inLocal[]){ + multipole = inMultipole; + local = inLocal; + } + + const FComplexe* getMultipole() const{ + return multipole; + } + + const FComplexe* getLocal() const{ + return local; + } + + FComplexe* getMultipole(){ + return multipole; + } + + FComplexe* getLocal(){ + return local; + } +}; + + +static const int P = 5; + +typedef FRotationParticle KernelParticleClass; +typedef FVector<KernelParticleClass> KernelContainerClass; +typedef KernelCell KernelCellClass; +typedef FRotationKernel<KernelParticleClass, KernelCellClass, KernelContainerClass, P > KernelClass; + +struct ScalFmmKernelHandle { + KernelClass** kernel; + int potentialDataSize; + int fieldDataSize; + int nbthread; +}; + +/******* Allocation : ******/ +int FmmKernel_init(void *fmmCore, void **fmmKernel){ + ScalFmmKernelHandle* kernelhandle = new ScalFmmKernelHandle; + memset(kernelhandle, 0, sizeof(ScalFmmKernelHandle)); + + int NbLevels; + FmmCore_getParameter(fmmCore, FMMCORE_TREE_HEIGHT, &NbLevels); + FReal boxWidth; + FmmCore_getParameter(fmmCore, FMMCORE_ROOT_BOX_WIDTH, &boxWidth); + FReal centerOfBox[3]; + FmmCore_getParameter(fmmCore, FMMCORE_ROOT_BOX_CENTER, centerOfBox); + + FmmCore_getParameter(fmmCore, FMMCORE_THREADS_NUMBER, &kernelhandle->nbthread); + + KernelClass original( NbLevels, boxWidth, FPoint(centerOfBox) ); + kernelhandle->kernel = new KernelClass*[kernelhandle->nbthread]; + for(int idxThread = 0 ; idxThread < kernelhandle->nbthread ; ++idxThread){ + kernelhandle->kernel[idxThread] = new KernelClass(original); + } + + kernelhandle->potentialDataSize = 1; + kernelhandle->fieldDataSize = 4; + + *fmmKernel = kernelhandle; + + return FMMAPI_NO_ERROR; +}/* : alloue et initialise le FmmKernel */ +int FmmKernel_free(void *fmmKernel){ + ScalFmmKernelHandle* kernelhandle = (ScalFmmKernelHandle*) fmmKernel; + + for(int idxThread = 0 ; idxThread < kernelhandle->nbthread ; ++idxThread){ + delete kernelhandle->kernel[idxThread]; + } + + delete[] kernelhandle->kernel; + delete kernelhandle; + + return FMMAPI_NO_ERROR; +} /* libére le FmmKernel */ + + + + +/******* Configuration : ***/ + +int FmmKernel_isParameterUsed(void * /*fmm*/, int *name, int *flag){ + switch( *name ){ + case FMMKERNEL_POTENTIAL_DATA_SIZE: + case FMMKERNEL_FIELD_DATA_SIZE: + case FMMKERNEL_HANDLES_P2P: + *flag = FMMAPI_SUPPORTED_PARAMETER; + break; + case FMMKERNEL_ACCURACY : + *flag = FMMAPI_UNSUPPORTED_PARAMETER; + break; + default: + *flag = FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_setParameter(void *fmmKernel, int *name, void*value){ + /*ScalFmmKernelHandle* kernelhandle = (ScalFmmKernelHandle*) fmmKernel;*/ + int flag; + + FmmKernel_isParameterUsed(fmmKernel, name, &flag); + if( flag != FMMAPI_SUPPORTED_PARAMETER){ + return flag; + } + + switch( *name ){ + case FMMKERNEL_POTENTIAL_DATA_SIZE : + case FMMKERNEL_FIELD_DATA_SIZE : + return FMMAPI_SUPPORTED_PARAMETER; + default: + return FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_setParameter(void *fmmKernel, int name, void*value){ + return FmmKernel_setParameter( fmmKernel, &name, value); +} + +int FmmKernel_getParameter(void *fmmKernel, int *name, void*value){ + ScalFmmKernelHandle* kernelhandle = (ScalFmmKernelHandle*) fmmKernel; + int flag; + + FmmKernel_isParameterUsed(fmmKernel, name, &flag); + if( flag != FMMAPI_SUPPORTED_PARAMETER){ + return flag; + } + + switch( *name ){ + case FMMKERNEL_POTENTIAL_DATA_SIZE : + *(int*)value = kernelhandle->potentialDataSize*sizeof(FReal); + break; + case FMMKERNEL_FIELD_DATA_SIZE : + *(int*)value = kernelhandle->fieldDataSize*sizeof(FReal); + break; + default: + return FMMAPI_UNKNOWN_PARAMETER; + } + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_getParameter(void *fmmKernel, int name, void*value){ + return FmmKernel_getParameter(fmmKernel, &name, value); +} + + +/****** Données FMM : *****/ +int FmmKernel_getMultipoleArraySize(void */*fmmCore*/, int *size) { + *size = ((P+2)*(P+1))/2 * sizeof(FComplexe); + return FMMAPI_NO_ERROR; +} /* Renvoie dans size la taille (en octets) de l'expansion multipôle associée à la boîte boxId */ + +int FmmKernel_getLocalArraySize(void */*fmmCore*/, int *size){ + *size = ((P+2)*(P+1))/2 * sizeof(FComplexe); + return FMMAPI_NO_ERROR; +} /* Renvoie dans size la taille (en octets) de l'expansion locale associée à la boîte boxId*/ + + +/******* Opérateurs FMM : **/ +int FmmKernel_P2M(void *fmmCore, void* boxId){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FComplexe* multipole; + FmmCore_getMultipoleArray(fmmCore, boxId, (void**)&multipole); + + KernelCellClass cell; + cell.attachArrays(multipole, 0); + int coord[3]; + FmmCore_getCoord(fmmCore, boxId, coord); + cell.setCoordinate(coord[0], coord[1], coord[2]); + + FReal* positions; + FReal* potentials; + int number; + FmmCore_getSource(fmmCore, boxId, &positions, (void**)&potentials, &number); + + KernelContainerClass sources; + KernelParticleClass part; + for(int idxPart = 0 ; idxPart < number ; ++idxPart){ + part.setPosition(FPoint(&positions[idxPart*3])); + part.setPhysicalValue(potentials[idxPart]); + sources.push(part); + } + + kernelhandle->kernel[threadId]->P2M(&cell, &sources); + + FmmCore_releaseSource(fmmCore, boxId, potentials, positions); + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_L2P(void *fmmCore, void* boxId){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FComplexe* local; + FmmCore_getLocalArray(fmmCore, boxId, (void**)&local); + + KernelCellClass cell; + cell.attachArrays(0,local); + int coord[3]; + FmmCore_getCoord(fmmCore, boxId, coord); + cell.setCoordinate(coord[0], coord[1], coord[2]); + + FReal* potentials; + FReal* positions; + int number; + //FmmCore_getTargetPoints(fmmCore, boxId, &positions, &number); + FmmCore_getSource(fmmCore, boxId, &positions, (void**)&potentials, &number); + + FReal* fields = new FReal[number*kernelhandle->fieldDataSize]; + FmmCore_getTargetField(fmmCore, boxId, fields); + + KernelContainerClass targets; + KernelParticleClass part; + for(int idxPart = 0 ; idxPart < number ; ++idxPart){ + part.setPosition(FPoint(&positions[idxPart*3])); + part.setForces(fields[idxPart*kernelhandle->fieldDataSize],fields[idxPart*kernelhandle->fieldDataSize+1], + fields[idxPart*kernelhandle->fieldDataSize+2]); + part.setPotential(fields[idxPart*kernelhandle->fieldDataSize+3]); + part.setPhysicalValue(potentials[idxPart]); + targets.push(part); + } + + kernelhandle->kernel[threadId]->L2P(&cell, &targets); + + for(int idxPart = 0 ; idxPart < number ; ++idxPart){ + fields[idxPart*kernelhandle->fieldDataSize] = targets[idxPart].getForces().getX(); + fields[idxPart*kernelhandle->fieldDataSize+1] = targets[idxPart].getForces().getY(); + fields[idxPart*kernelhandle->fieldDataSize+2] = targets[idxPart].getForces().getZ(); + fields[idxPart*kernelhandle->fieldDataSize+3] = targets[idxPart].getPotential(); + } + + FmmCore_releaseTargetPoints(fmmCore, boxId, positions); + FmmCore_setTargetField(fmmCore, boxId, fields); + delete[] fields; + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_M2M(void *fmmCore, void *boxIdFather, void *boxIdSon){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FComplexe* multipole; + FmmCore_getMultipoleArray(fmmCore, boxIdFather, (void**)&multipole); + + KernelCellClass cellFather; + cellFather.attachArrays(multipole, 0); + int coordFather[3]; + FmmCore_getCoord(fmmCore, boxIdFather, coordFather); + cellFather.setCoordinate(coordFather[0], coordFather[1], coordFather[2]); + + FmmCore_getMultipoleArray(fmmCore, boxIdSon, (void**)&multipole); + + KernelCellClass cellSon; + cellSon.attachArrays(multipole, 0); + int coordChild[3]; + FmmCore_getCoord(fmmCore, boxIdSon, coordChild); + cellSon.setCoordinate(coordChild[0], coordChild[1], coordChild[2]); + + int level; + FmmCore_getLevel(fmmCore,boxIdFather, &level); + + const KernelCellClass* children[8]; + memset(children, 0, sizeof(KernelCellClass*)*8); + const int mindex = ((coordChild[0]&1) * 2 + (coordChild[1]&1)) * 2 + (coordChild[2]&1); + children[mindex] = &cellSon; + + kernelhandle->kernel[threadId]->M2M(&cellFather, children, level); + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_L2L(void *fmmCore, void *boxIdFather, void *boxIdSon){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FComplexe* local; + FmmCore_getLocalArray(fmmCore, boxIdFather, (void**)&local); + + KernelCellClass cellFather; + cellFather.attachArrays(0, local); + int coordFather[3]; + FmmCore_getCoord(fmmCore, boxIdFather, coordFather); + cellFather.setCoordinate(coordFather[0], coordFather[1], coordFather[2]); + + FmmCore_getLocalArray(fmmCore, boxIdSon, (void**)&local); + + KernelCellClass cellSon; + cellSon.attachArrays(0, local); + int coordChild[3]; + FmmCore_getCoord(fmmCore, boxIdSon, coordChild); + cellSon.setCoordinate(coordChild[0], coordChild[1], coordChild[2]); + + int level; + FmmCore_getLevel(fmmCore,boxIdFather, &level); + + KernelCellClass* children[8]; + memset(children, 0, sizeof(KernelCellClass*)*8); + const int mindex = ((coordChild[0]&1) * 2 + (coordChild[1]&1)) * 2 + (coordChild[2]&1); + children[mindex] = &cellSon; + + kernelhandle->kernel[threadId]->L2L(&cellFather, children, level); + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_M2L(void *fmmCore, void *boxIdSrc, void *boxIdDest){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FComplexe* multipole; + FmmCore_getMultipoleArray(fmmCore, boxIdSrc, (void**)&multipole); + KernelCellClass cellSrc; + cellSrc.attachArrays(multipole,0); + int coord[3]; + FmmCore_getCoord(fmmCore, boxIdSrc, coord); + cellSrc.setCoordinate(coord[0], coord[1], coord[2]); + + FComplexe* local; + FmmCore_getLocalArray(fmmCore, boxIdDest, (void**)&local); + KernelCellClass cellDst; + cellDst.attachArrays(0, local); + FmmCore_getCoord(fmmCore, boxIdDest, coord); + cellDst.setCoordinate(coord[0], coord[1], coord[2]); + + int level; + FmmCore_getLevel(fmmCore, boxIdDest, &level); + + const int xdiff = cellSrc.getCoordinate().getX() - cellDst.getCoordinate().getX(); + const int ydiff = cellSrc.getCoordinate().getY() - cellDst.getCoordinate().getY(); + const int zdiff = cellSrc.getCoordinate().getZ() - cellDst.getCoordinate().getZ(); + const int index = (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3; + + const KernelCellClass* inter[343]; + memset(inter, 0, sizeof(KernelCellClass*)*343); + inter[index] = &cellSrc; + + kernelhandle->kernel[threadId]->M2L(&cellDst, inter, 1, level); + + return FMMAPI_NO_ERROR; +} + +int FmmKernel_P2P(void *fmmCore, void *boxIdSrc, void *boxIdDest){ + ScalFmmCoreHandle* corehandle = (ScalFmmCoreHandle*)fmmCore; + ScalFmmKernelHandle* kernelhandle; + FmmCore_getKernelData(corehandle, (void**)&kernelhandle); + int threadId; + FmmCore_getParameter(fmmCore, FMMCORE_THREAD_ID, &threadId); + + FReal* positionsTargets; + FReal* potentialsTargets; + int numberTargets; + //FmmCore_getTargetPoints(fmmCore, boxIdDest, &positionsTargets, &numberTargets); + FmmCore_getSource(fmmCore, boxIdDest, &positionsTargets, (void**)&potentialsTargets, &numberTargets); + + FReal* fieldsTargets = new FReal[numberTargets*kernelhandle->fieldDataSize]; + FmmCore_getTargetField(fmmCore, boxIdDest, fieldsTargets); + + int coordTargets[3]; + FmmCore_getCoord(fmmCore, boxIdDest, coordTargets); + + KernelContainerClass targets; + KernelParticleClass part; + for(int idxPart = 0 ; idxPart < numberTargets ; ++idxPart){ + part.setPosition(FPoint(&positionsTargets[idxPart*3])); + part.setForces(fieldsTargets[idxPart*kernelhandle->fieldDataSize], + fieldsTargets[idxPart*kernelhandle->fieldDataSize+1], + fieldsTargets[idxPart*kernelhandle->fieldDataSize+2]); + part.setPotential(fieldsTargets[idxPart*kernelhandle->fieldDataSize+3]); + part.setPhysicalValue(potentialsTargets[idxPart]); + targets.push(part); + } + + FReal* positionsSrc; + FReal* potentialsSrc; + int numberSources; + FmmCore_getSource(fmmCore, boxIdSrc, &positionsSrc, (void**)&potentialsSrc, &numberSources); + + KernelContainerClass sources; + for(int idxPart = 0 ; idxPart < numberSources ; ++idxPart){ + part.setPosition(FPoint(&positionsSrc[idxPart*3])); + part.setPhysicalValue(potentialsSrc[idxPart]); + sources.push(part); + } + + //KernelContainerClass* neigh[27]; + //memset(neigh, 0, sizeof(KernelContainerClass*)*27); + //neigh[0] = &sources; + //kernelhandle->kernel->P2PRemote(FTreeCoordinate(coordTargets[0],coordTargets[1],coordTargets[2]), + // &targets, &sources, neigh, 1); + for(int idxTarget = 0 ; idxTarget < numberTargets ; ++idxTarget){ + for(int idxSource = 0 ; idxSource < numberSources ; ++idxSource){ + FReal dx = sources[idxSource].getPosition().getX() - targets[idxTarget].getPosition().getX(); + FReal dy = sources[idxSource].getPosition().getY() - targets[idxTarget].getPosition().getY(); + FReal dz = sources[idxSource].getPosition().getZ() - targets[idxTarget].getPosition().getZ(); + const FReal distSquare = (dx*dx + dy*dy + dz*dz); + + if(distSquare > 10E-6*10E-6){ + kernelhandle->kernel[threadId]->particlesInteraction(&targets[idxTarget],sources[idxSource]); + } + } + } + + for(int idxPart = 0 ; idxPart < numberTargets ; ++idxPart){ + fieldsTargets[idxPart*kernelhandle->fieldDataSize] = targets[idxPart].getForces().getX(); + fieldsTargets[idxPart*kernelhandle->fieldDataSize+1] = targets[idxPart].getForces().getY(); + fieldsTargets[idxPart*kernelhandle->fieldDataSize+2] = targets[idxPart].getForces().getZ(); + fieldsTargets[idxPart*kernelhandle->fieldDataSize+3] = targets[idxPart].getPotential(); + } + + FmmCore_releaseSource(fmmCore, boxIdDest, potentialsSrc, positionsSrc); + FmmCore_releaseTargetPoints(fmmCore, boxIdDest, positionsTargets); + FmmCore_setTargetField(fmmCore, boxIdDest, fieldsTargets); + delete[] fieldsTargets; + + return FMMAPI_NO_ERROR; +} /* pas mutuel, i.e. on fait seulement dans 1 sens. */ + + diff --git a/Addons/FmmApi/Tests/testFmmApi.cpp b/Addons/FmmApi/Tests/testFmmApi.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0aa63db719e8454547f87b8306c61af8bbb7843b --- /dev/null +++ b/Addons/FmmApi/Tests/testFmmApi.cpp @@ -0,0 +1,192 @@ +// =================================================================================== +// 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 <iostream> +#include <cstdio> + +#include "../../Src/Utils/FParameters.hpp" +#include "../../Src/Utils/FTic.hpp" + +#include "../../Src/Files/FRandomLoader.hpp" + + +/** This program show an example of use of the fmm api + */ + +// Simply create particles and try the kernels +int main(int argc, char ** argv){ + ///////////////////////What we do///////////////////////////// + std::cout << ">> This executable has to be used to test the FMM API\n"; + ////////////////////////////////////////////////////////////// + + int NbLevels = FParameters::getValue(argc,argv,"-h", 7); + int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3); + const int NbPart = FParameters::getValue(argc,argv,"-nb", 2000000); + FTic counter; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + FRandomLoader<KernelParticleClass> loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1); + + void* FmmCoreHandle; + FmmCore_init(&FmmCoreHandle); + + FmmCore_setParameter(FmmCoreHandle, FMMCORE_TREE_HEIGHT, &NbLevels); + FReal boxWidth = loader.getBoxWidth(); + FmmCore_setParameter(FmmCoreHandle, FMMCORE_ROOT_BOX_WIDTH, &boxWidth); + FmmCore_setParameter(FmmCoreHandle, FMMCORE_ROOT_BOX_CENTER, loader.getCenterOfBox().getDataValue()); + + FmmCore_finishInit(FmmCoreHandle); + + void* FmmKernelHandle; + FmmKernel_init(FmmCoreHandle, &FmmKernelHandle); + + FmmCore_setKernelData(FmmCoreHandle, FmmKernelHandle); + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl; + std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl; + counter.tic(); + + int nbPart = (int)loader.getNumberOfParticles(); + FReal* potentials = new FReal[nbPart]; + FReal* positions = new FReal[nbPart*3]; + + { + KernelParticleClass part; + part.setPhysicalValue(0.1); + + for(int idx = 0 ; idx < nbPart ; ++idx){ + loader.fillParticle(part); + + potentials[idx] = part.getPhysicalValue(); + positions[3*idx] = part.getPosition().getX(); + positions[3*idx+1] = part.getPosition().getY(); + positions[3*idx+2] = part.getPosition().getZ(); + } + + FmmCore_setPositions(FmmCoreHandle, &nbPart, positions); + FmmCore_setPotentials(FmmCoreHandle, potentials); + } + + counter.tac(); + std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Working on particles ..." << std::endl; + counter.tic(); + + FmmCore_doComputation(FmmCoreHandle); + + counter.tac(); + std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Getting results ..." << std::endl; + counter.tic(); + + FReal* fields = new FReal[4*nbPart]; + + FmmCore_getField(FmmCoreHandle, fields); + + counter.tac(); + std::cout << "Done " << "(@retrieve = " << counter.elapsed() << "s)." << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Computing direct ..." << std::endl; + counter.tic(); + + FReal* fieldsdirect = new FReal[4*nbPart]; + FMemUtils::setall(fieldsdirect, 0.0, nbPart * 4); + + for(int idxTarget = 0 ; idxTarget < nbPart ; ++idxTarget){ + + for(int idxSource = idxTarget + 1 ; idxSource < nbPart ; ++idxSource){ + FReal dx = positions[idxSource*3] - positions[idxTarget*3]; + FReal dy = positions[idxSource*3+1] - positions[idxTarget*3+1]; + FReal dz = positions[idxSource*3+2] - positions[idxTarget*3+2]; + + FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz); + FReal inv_distance = FMath::Sqrt(inv_square_distance); + + inv_square_distance *= inv_distance; + inv_square_distance *= potentials[idxTarget] * potentials[idxSource]; + + dx *= inv_square_distance; + dy *= inv_square_distance; + dz *= inv_square_distance; + + fieldsdirect[4*idxTarget] += dx; + fieldsdirect[4*idxTarget+1] += dy; + fieldsdirect[4*idxTarget+2] += dz; + fieldsdirect[4*idxTarget+3] += inv_distance * potentials[idxSource]; + + fieldsdirect[4*idxSource] += -dx; + fieldsdirect[4*idxSource+1] += -dy; + fieldsdirect[4*idxSource+2] += -dz; + fieldsdirect[4*idxSource+3] += inv_distance * potentials[idxTarget]; + } + } + + + counter.tac(); + std::cout << "Done " << "(@direct = " << counter.elapsed() << "s)." << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + + std::cout << "Comparing results ..." << std::endl; + counter.tic(); + + FMath::FAccurater forces; + FMath::FAccurater potential; + + for(int idx = 0 ; idx < nbPart ; ++idx){ + if(idx < 10){ + std::cout << fieldsdirect[4*idx] << " fmm " << fields[4*idx] << std::endl; + } + forces.add( fieldsdirect[4*idx] ,fields[4*idx]); + forces.add( fieldsdirect[4*idx+1] ,fields[4*idx+1]); + forces.add( fieldsdirect[4*idx+2] ,fields[4*idx+2]); + potential.add( fieldsdirect[4*idx+3] ,fields[4*idx+3]); + } + + counter.tac(); + std::cout << "Done " << "(@comparing = " << counter.elapsed() << "s)." << std::endl; + + std::cout << "\tForces inf " << forces.getInfNorm() << " normL2 " << forces.getL2Norm() << std::endl; + std::cout << "\tPotential inf " << potential.getInfNorm() << " normL2 " << potential.getL2Norm() << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + + delete[] fieldsdirect; + delete[] fields; + delete[] potentials; + delete[] positions; + + FmmCore_free(FmmCoreHandle); + FmmKernel_free(FmmKernelHandle); + + return 0; +} + + +