Commit aa4c4d19 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

API can now call chebyshev/lagrange kernel succesfully (only cheb has been...

API can now call chebyshev/lagrange kernel succesfully (only cheb has been tested), remaining to do : updates positions of particles
parent db0c2164
......@@ -68,6 +68,29 @@ typedef enum kernel_type {
} scalfmm_kernel_type;
/**
* Enum over the different way scalfmm can handle a particule moving
* out of the simulation box
*/
typedef enum out_of_the_box_config {
exiting = 0, /*If a particule move outside the simulation box,
the simulation stop*/
periodic = 1, /*If a particule move outside the simulation box,
the particules is inserted back at the other side
of the simulation box*/
erasing = 2 /*If a particule move outside the simulation box, it
simply disappears from the simulation */
} scalfmm_out_of_box_behavior;
/**
* Enum over the different algorithm that scalfmm can use
*/
typedef enum scalfmm_algorithm_config {
sequential = 0, /* Use the sequential version of Scalfmm*/
multi_thread = 1 /* Use the Multi thread version of Scalfmm*/
} scalfmm_algorithm;
/**
* Handle for the user
*/
......@@ -135,6 +158,11 @@ void scalfmm_tree_insert_particles_xyz(scalfmm_handle Handle, int NbPositions, d
* the array. First particle inserted will take value physicalValues[0].
*/
void scalfmm_set_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues);
/**
* @brief get the physical values.
*
* WARNING : the user must allocate (and initialize) the array given
*/
void scalfmm_get_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues);
/**
......@@ -206,7 +234,7 @@ void scalfmm_set_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfPart
* array.
*/
void scalfmm_get_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToFill);
void scalfmm_set_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToFill);
void scalfmm_set_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToRead);
void scalfmm_get_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToFill);
void scalfmm_set_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToFill);
......@@ -251,7 +279,24 @@ void scalfmm_get_positions(scalfmm_handle Handle, int NbPositions, double * X, d
void scalfmm_get_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * updatedXYZ);
void scalfmm_get_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z);
/**
* @brief This function provides a way for the user to define scalfmm
* behavior in case a particule get out of the box after a
* displacement
* @param Handle scalfmm_handle provided by scalfmm_init
* @param Member of enum scalfmm_out_of_box_behavior
*/
void scalfmm_out_of_the_box_config(scalfmm_handle Handle,scalfmm_out_of_box_behavior config);
/**
* @brief This function provides a way for choosing the algorithm to
* be used
* @param Handle scalfmm_handle provided by scalfmm_init
* @param Member of enum scalfmm_algorithm
*/
void scalfmm_algorithm_config(scalfmm_handle Handle,scalfmm_algorithm config);
/////////////////////////////////////////////////////////////////////
......
......@@ -25,11 +25,10 @@
#include "FScalFMMEngine.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
/**
......@@ -38,31 +37,363 @@
*/
template<class InterCell,class InterKernel,
class ContainerClass = FP2PParticleContainerIndexed<>,
class LeafClass = FSimpleLeaf<FP2PParticleContainerIndexed<>>,
class LeafClass = FSimpleLeaf<FP2PParticleContainerIndexed<> >,
class MatrixKernelClass = FInterpMatrixKernelR>
class FInterEngine : public FScalFMMEngine{
private:
//Pointer to the kernel to be executed
InterKernel * kernel;
MatrixKernelClass * matrix;
//Link to the tree
FOctree<InterCell,ContainerClass,LeafClass> * octree;
FOctree<InterCell,ContainerClass,LeafClass> * octree;
public:
/**
* @brief Constructor : its only build the tree
* @brief Constructor : build the tree and the interpolation
* kernel
* @param TreeHeight Height of the tree
* @param BoxWidth box Width
* @param BoxCenter double[3] coordinate of the center of the
* simulation box
*/
FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter) : kernel(nullptr), octree(nullptr){
FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter) : kernel(nullptr), matrix(nullptr), octree(nullptr){
octree = new FOctree<InterCell,ContainerClass,LeafClass>(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint(BoxCenter));
this->matrix = new MatrixKernelClass();
this->kernel = new InterKernel(TreeHeight,BoxWidth,FPoint(BoxCenter),matrix);
kernelType = chebyshev;
}
//TODO free kernel too
~FInterEngine(){
free(octree);
free(kernel);
}
//Inserting array of position
void tree_insert_particles_xyz(int NbPositions, double * XYZ){
for(int idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint(&XYZ[3*idPart]),idPart);
}
nbPart += NbPositions;
}
//Inserting arrayS of position
void tree_insert_particles(int NbPositions, double * X, double * Y, double * Z){
for(int idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint(X[idPart],Y[idPart],Z[idPart]),idPart);
}
nbPart += NbPositions;
}
//Set the physical values
void set_physical_values(int nbPhysicalValues,double * physicalValues){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
sources->getPhysicalValues()[idxPart] = physicalValues[indexes[idxPart]];
}
});
}
//Set only a subpart of physical values
//Algorithm : loop over each leaf, and then search in user array
//if any index matches
void set_physical_values_npart( int nbPhysicalValues, int* idxOfParticles, double * physicalValues){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbPhysicalValues && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
sources->getPhysicalValues()[idxPart] = physicalValues[indexes[idxPart]];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
//get back the physical values
void get_physical_values( int nbPhysicalValues, double * physicalValues){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
physicalValues[indexes[idxPart]] = sources->getPhysicalValues()[idxPart];
}
});
}
//Same algorithm as in set_physical_values_npart
void get_physical_values_npart( int nbPhysicalValues, int* idxOfParticles, double * physicalValues){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbPhysicalValues && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
physicalValues[indexes[idxPart]] = sources->getPhysicalValues()[idxPart];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
void get_forces_xyz( int nbParts, double * forcesToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
}
});
}
void get_forces_xyz_npart(int nbParts, int* idxOfParticles , double * forcesToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbParts && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
void get_forces( int nbParts, double * fX, double* fY, double* fZ){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
fX[indexes[idxPart]] = sources->getForcesX()[idxPart];
fY[indexes[idxPart]] = sources->getForcesY()[idxPart];
fZ[indexes[idxPart]] = sources->getForcesZ()[idxPart];
}
});
}
void get_forces_npart(int nbParts, int* idxOfParticles ,double * fX, double* fY, double* fZ){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbParts && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
fX[indexes[idxPart]] = sources->getForcesX()[idxPart];
fY[indexes[idxPart]] = sources->getForcesY()[idxPart];
fZ[indexes[idxPart]] = sources->getForcesZ()[idxPart];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
//To set initial condition
void set_forces_xyz( int nbParts, double * forcesToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
sources->getForcesX()[idxPart] = forcesToFill[indexes[idxPart]*3+0];
sources->getForcesY()[idxPart] = forcesToFill[indexes[idxPart]*3+1];
sources->getForcesZ()[idxPart] = forcesToFill[indexes[idxPart]*3+2];
}
});
}
void set_forces_xyz_npart( int nbParts, int* idxOfParticles, double * forcesToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbParts && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
sources->getForcesX()[idxPart] = forcesToFill[indexes[idxPart]*3+0];
sources->getForcesY()[idxPart] = forcesToFill[indexes[idxPart]*3+1];
sources->getForcesZ()[idxPart] = forcesToFill[indexes[idxPart]*3+2];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
void set_forces( int nbParts, double * fX, double* fY, double* fZ){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
sources->getForcesX()[idxPart] = fX[indexes[idxPart]];
sources->getForcesY()[idxPart] = fY[indexes[idxPart]];
sources->getForcesZ()[idxPart] = fZ[indexes[idxPart]];
}
});
}
void set_forces_npart( int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbParts && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
sources->getForcesX()[idxPart] = fX[indexes[idxPart]];
sources->getForcesY()[idxPart] = fY[indexes[idxPart]];
sources->getForcesZ()[idxPart] = fZ[indexes[idxPart]];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
//Set the potentials
void set_potentials(int nbPotentials,double * potentialsToRead){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
sources->getPotentials()[idxPart] = potentialsToRead[indexes[idxPart]];
}
});
}
//Set only a subpart of potentials
//Algorithm : loop over each leaf, and then search in user array
//if any index matches
void set_potentials_npart( int nbPotentials, int* idxOfParticles, double * potentialsToRead){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbPotentials && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
sources->getPotentials()[idxPart] = potentialsToRead[indexes[idxPart]];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
//get back the potentials
void get_potentials( int nbPotentials, double * potentialsToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
potentialsToFill[indexes[idxPart]] = sources->getPotentials()[idxPart];
}
});
}
//Same algorithm as in set_potentials_npart
void get_potentials_npart( int nbPotentials, int* idxOfParticles, double * potentialsToFill){
octree->forEachLeaf([&](LeafClass* leaf){
ContainerClass * sources = leaf->getSrc();
FVector<int> indexes = sources->getIndexes();
int nbPartThere = sources->getNbParticles();
for(int idxPart = 0 ; idxPart<nbPartThere ; ++idxPart){
int iterPart = 0;
bool notFoundYet = true;
while(iterPart < nbPotentials && notFoundYet){
if(indexes[idxPart] == idxOfParticles[iterPart]){
potentialsToFill[indexes[idxPart]] = sources->getPotentials()[idxPart];
notFoundYet = false;
}
else{
++iterPart;
}
}
}
});
}
// class ContainerClass = FP2PParticleContainerIndexed<>,
// class LeafClass = FSimpleLeaf<FP2PParticleContainerIndexed<> >,
// class MatrixKernelClass = FInterpMatrixKernelR>
void execute_fmm(){
typedef FOctree<InterCell,ContainerClass,LeafClass> OctreeClass;
switch(Algorithm){
case 0:
{
typedef FFmmAlgorithm<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassSeq;
AlgoClassSeq algoSeq(octree,kernel);
algoSeq.execute();
break;
}
case 1:
{
typedef FFmmAlgorithmThread<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassThread;
AlgoClassThread algoThread(octree,kernel);
algoThread.execute();
break;
}
default :
std::cout<< "No algorithm found (probably for strange reasons) : "<< Algorithm <<" exiting" << std::endl;
}
}
};
......
......@@ -47,17 +47,39 @@
* @class FScalFMMEngine
*/
class FScalFMMEngine{
private:
scalfmm_kernel_type kernel;
protected:
scalfmm_kernel_type kernelType;
scalfmm_out_of_box_behavior OutOfBoxBehavior;
scalfmm_algorithm Algorithm;
FVector<bool> progress;
int nbPart;
public:
FScalFMMEngine() : OutOfBoxBehavior(exiting), 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
//kernel
scalfmm_kernel_type getKernelType(){
return this->kernel;
return this->kernelType;
}
//To deal with particles moving outside the box
void out_of_the_box_config(scalfmm_out_of_box_behavior config){
this->OutOfBoxBehavior = config;
}
//Function about the tree (to be redefined)
//To change default algorithm
void algorithm_config(scalfmm_algorithm config){
this->Algorithm = config;
}
//Functions displayed there are function that are to be redefined
//by specific Engine
//Function about the tree
virtual void tree_insert_particles( int NbPositions, double * arrayX, double * arrayY, double * arrayZ){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
......@@ -96,6 +118,7 @@ public:
virtual void get_forces_npart( int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
//To set initial condition
virtual void set_forces_xyz( int nbParts, double * forcesToFill){
FAssertLF(0,"No tree instancied, exiting ...\n");
......@@ -113,17 +136,17 @@ public:
virtual void get_potentials( int nbParts, double * potentialsToFill){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void set_potentials( int nbParts, double * potentialsToFill){
virtual void set_potentials( int nbParts, double * potentialsToRead){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void get_potentials_npart( int nbParts, int* idxOfParticles, double * potentialsToFill){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void set_potentials_npart( int nbParts, int* idxOfParticles, double * potentialsToFill){
virtual void set_potentials_npart( int nbParts, int* idxOfParticles, double * potentialsToRead){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
//To deal with positions
//Function to move particles
virtual void add_to_positions_xyz( int NbPositions, double * updatedXYZ){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
......@@ -134,7 +157,7 @@ public:
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void add_to_positions_npart( int NbPositions, int* idxOfParticles,
double * X, double * Y , double * Z){
double * X, double * Y , double * Z){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void set_positions_xyz( int NbPositions, double * updatedXYZ){
......@@ -147,7 +170,7 @@ public:
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void set_positions_npart( int NbPositions, int* idxOfParticles,
double * X, double * Y , double * Z){
double * X, double * Y , double * Z){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void get_positions_xyz( int NbPositions, double * updatedXYZ){
......@@ -160,7 +183,7 @@ public:
FAssertLF(0,"No tree instancied, exiting ...\n");
}
virtual void get_positions_npart( int NbPositions, int* idxOfParticles,
double * X, double * Y , double * Z){
double * X, double * Y , double * Z){
FAssertLF(0,"No tree instancied, exiting ...\n");
}
......@@ -170,8 +193,8 @@ public:
FAssertLF(0,"No user kernel defined, exiting ...\n");
}
virtual void scalfmm_execute_fmm(){
FAssertLF(0,"No kernel set, exiting ...\n");
virtual void execute_fmm(){
FAssertLF(0,"No kernel set, cannot execute anything, exiting ...\n");
}
};
......@@ -272,6 +295,11 @@ extern "C" void scalfmm_set_potentials_npart(scalfmm_handle Handle, int nbParts,
//To deal with positions
//Out of the box behavior
extern "C" void scalfmm_out_of_the_box_config(scalfmm_handle Handle,scalfmm_out_of_box_behavior config){
((ScalFmmCoreHandle * ) Handle)->engine->out_of_the_box_config(config);
}
//Update
extern "C" void scalfmm_add_to_positions_xyz(scalfmm_handle Handle, int NbPositions, double * updatedXYZ){
((ScalFmmCoreHandle * ) Handle)->engine->add_to_positions_xyz(NbPositions, updatedXYZ);
......@@ -286,7 +314,7 @@ extern "C" void scalfmm_add_to_positions_xyz_npart(scalfmm_handle Handle, int Nb
}
extern "C" void scalfmm_add_to_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles,
double * X, double * Y , double * Z){
double * X, double * Y , double * Z){
((ScalFmmCoreHandle * ) Handle)->engine->add_to_positions_npart(NbPositions, idxOfParticles, X, Y, Z);
}
//Set new positions
......@@ -320,14 +348,20 @@ extern "C" void scalfmm_get_positions_xyz_npart(scalfmm_handle Handle, int NbPos
}
extern "C" void scalfmm_get_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles,
double * X, double * Y , double * Z){
double * X, double * Y , double * Z){
((ScalFmmCoreHandle * ) Handle)->engine->get_positions_npart(NbPositions, idxOfParticles, X, Y, Z);
}
//To choose algorithm
extern "C" void scalfmm_algorithm_config(scalfmm_handle Handle, scalfmm_algorithm config){
((ScalFmmCoreHandle * ) Handle)->engine->algorithm_config(config);
}