Commit b8eec341 authored by COULAUD Olivier's avatar COULAUD Olivier

Add periodic boundary conditions in Exemples ; Check periodic boundary conditions.

parent e089d8b5
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// Keep in private GIT
//
//
/** \brief Chebyshev FMM example
*
* \file
* \authors O. Coulaud
*
* This program runs the FMM Algorithm with the interpolation kernel based on
* Chebyshev (grid points) interpolation (1/r kernel). It then compares the
* results with a direct computation.
*/
#include <string>
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
//
template<typename FReal, int ORDER>
using FInterpolationCell = FChebCell<FReal, ORDER>;
template<typename FReal, typename GroupCellClass,
typename GroupContainerClass,
typename MatrixKernelClass, int ORDER>
using FInterpolationKernel = FChebSymKernel<FReal,
GroupCellClass,
GroupContainerClass,
MatrixKernelClass,
ORDER> ;
const std::string interpolationType("Chebyshev interpolation");
#include "MPIInterpolationFMM.hpp"
// See LICENCE file at project root
// ==== CMAKE =====
// @FUSE_MPI
// @FUSE_BLAS
// ================
#include <iostream>
#include <stdexcept>
#include <cstdio>
#include <cstdlib>
#include "ScalFmmConfig.h"
#include "Containers/FOctree.hpp"
#include "Utils/FMpi.hpp"
#include "Core/FFmmAlgorithmThreadProc.hpp"
#include "Core/FFmmAlgorithmThreadProcPeriodic.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
#include "Utils/FLeafBalance.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
//
// Order of the Interpolation approximation
static constexpr unsigned ORDER = 6 ;
using FReal = double;
// 1/r kernel
//
using MatrixKernelClass = FInterpMatrixKernelR<FReal> ;
//
// Typedefs
using ContainerClass = FP2PParticleContainerIndexed<FReal>;
using LeafClass = FSimpleLeaf<FReal, ContainerClass>;
using CellClass = FInterpolationCell<FReal, ORDER>;
using OctreeClass = FOctree<FReal,CellClass,ContainerClass,LeafClass>;
using MatrixKernelClass = FInterpMatrixKernelR<FReal>;
const MatrixKernelClass MatrixKernel;
using KernelClass = FInterpolationKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> ;
using FmmClassProc = FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
using FmmClassProcPER = FFmmAlgorithmThreadProcPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
/// \file
//!
//! \brief This program runs the MPI FMM with Chebyshev/Lagrange interpolation of 1/r kernel
//! \authors B. Bramas, O. Coulaud
//!
//! This code is a short example to use the FMM Algorithm Proc with Chebyshev or equispaced grid points Interpolation for the 1/r kernel
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
///////// PARAMETERS HANDLING //////////////////////////////////////
// const FParameterNames LocalOptionPeriodicDeep { {"-periodic","-per"}, "Perdiodic boundary condition (-per 5) the box grow by a facor (3L)^5"};
FHelpDescribeAndExit(argc, argv,
"Driver for Chebyshev Interpolation kernel using MPI (1/r kernel). "
"Usully run using : mpirun -np nb_proc_needed ./ChebyshevInterpolationAlgorithm [params].",
FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight,
FParameterDefinitions::InputFile,
FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads,
FParameterDefinitions::PeriodicityNbLevels);
const std::string defaultFile("../Data/test20k.fma");
const std::string filename = FParameters::getStr(argc,argv, FParameterDefinitions::InputFile.options, defaultFile.c_str());
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
bool periodicCondition = false ;
if(FParameters::existParameter(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options)){
periodicCondition = true;
}
const unsigned int aboveTree = FParameters::getValue(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options, 5);
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
//
std::cout << "Parameters"<< std::endl
<< " Octree Depth " << TreeHeight << std::endl
<< " SubOctree depth " << SubTreeHeight << std::endl;
if(periodicCondition){
std::cout << " AboveTree "<< aboveTree <<std::endl;
}
std::cout << " Input file name: " << filename << std::endl
<< " Thread count : " << NbThreads << std::endl
<< std::endl;
///////// VAR INIT /////////////////////////////////////////////////
// Initialize values for MPI
FMpi app(argc,argv);
//
// Initialize timer
FTic time;
// Creation of the particle loader
FMpiFmaGenericLoader<FReal> loader(filename,app.global());
if(!loader.isOpen()) {
throw std::runtime_error("Particle file couldn't be opened!") ;
}
// Initialize empty oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FSize localParticlesNumber = 0 ;
{ // -----------------------------------------------------
if(app.global().processId() == 0){
std::cout << "Loading & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
}
time.tic();
/* Mock particle structure to balance the tree over the processes. */
struct TestParticle{
FSize index; // Index of the particle in the original file.
FPoint<FReal> position; // Spatial position of the particle.
FReal physicalValue; // Physical value of the particle.
/* Returns the particle position. */
const FPoint<FReal>& getPosition(){
return position;
}
};
// Temporary array of particles read by this process.
TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
memset(particles, 0, (sizeof(TestParticle) * loader.getMyNumberOfParticles()));
// Index (in file) of the first particle that will be read by this process.
FSize idxStart = loader.getStart();
std::cout << "Proc:" << app.global().processId() << " start-index: " << idxStart << std::endl;
// Read particles from parts.
for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){
// Store the index (in the original file) the particle.
particles[idxPart].index = idxPart + idxStart;
// Read particle from file
loader.fillParticle(&particles[idxPart].position,
&particles[idxPart].physicalValue);
}
// Final vector of particles
FVector<TestParticle> finalParticles;
FLeafBalance balancer;
// Redistribute particules between processes
FMpiTreeBuilder< FReal, TestParticle >::
DistributeArrayToContainer(app.global(),
particles,
loader.getMyNumberOfParticles(),
tree.getBoxCenter(),
tree.getBoxWidth(),
tree.getHeight(),
&finalParticles,
&balancer);
// Free temporary array memory.
delete[] particles;
// Insert final particles into tree.
for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position,
finalParticles[idx].index,
finalParticles[idx].physicalValue);
}
time.tac();
localParticlesNumber = finalParticles.getSize() ;
double timeUsed = time.elapsed();
double minTime,maxTime;
std::cout << "Proc:" << app.global().processId()
<< " " << finalParticles.getSize()
<< " particles have been inserted in the tree. (@Reading and Inserting Particles = "
<< time.elapsed() << " s)."
<< std::endl;
MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
if(app.global().processId() == 0){
std::cout << "readinsert-time-min:" << minTime
<< " readinsert-time-max:" << maxTime
<< std::endl;
}
} // -----------------------------------------------------
std::vector<MortonIndex> mortonLeafDistribution(2*app.global().processCount());
{ // -----------------------------------------------------
std::cout << "\n"<<interpolationType<<" FMM Proc (ORDER="<< ORDER << ") ... " << std::endl;
time.tic();
// Kernels to use (pointer because of the limited size of the stack)
FAbstractAlgorithm * algorithm = nullptr;
FAlgorithmTimers * timer = nullptr;
// non periodic FMM algorithm
std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, loader.getBoxWidth(),
loader.getCenterOfBox(),&MatrixKernel));
FmmClassProc algoNoPer(app.global(),&tree, kernelsNoPer.get());
//
// periodic FMM algorithm
FmmClassProcPER algoPer(app.global(),&tree, aboveTree);
KernelClass kernelsPer(algoPer.extendedTreeHeight(), algoPer.extendedBoxWidth(),
algoPer.extendedBoxCenter(),&MatrixKernel);
algoPer.setKernel(&kernelsPer);
///////////////////////////////////////////////////////////////////////////////////////////////////
if(! periodicCondition) {// Non periodic case
algorithm = &algoNoPer ;
timer = &algoNoPer ;
}
else { // Periodic case
algorithm = &algoPer ;
timer = &algoPer ;
}
//
// FMM exectution FFmmFarField
algorithm->execute();
//
time.tac();
double timeUsed = time.elapsed();
double minTime,maxTime;
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl;
MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
if(app.global().processId() == 0){
std::cout << "exec-time-min:" << minTime
<< " exec-time-max:" << maxTime
<< std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ;
FReal energy =0.0 ;
//
// Loop over all leaves
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<" potential " << potentials[idxPart]
<< " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
});
FReal gloEnergy = app.global().reduceSum(energy);
if(0 == app.global().processId()){
std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
}
// -----------------------------------------------------
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
algorithm->getAndBuildMortonIndexAtLeaf(mortonLeafDistribution);
//
if(app.global().processId() == 0)
{
std::cout << " Morton distribution "<< std::endl ;
for(auto v : mortonLeafDistribution)
std::cout << v << " ";
std::cout << std::endl;
}
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FMpiFmaGenericWriter<FReal> paraWriter(name,app);
paraWriter.writeDistributionOfParticlesFromOctree(tree,loader.getNumberOfParticles(),localParticlesNumber,mortonLeafDistribution);
}
return 0;
}
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ================
// Keep in private GIT
//
//
/** \brief Uniform FMM example
*
* \file
* \authors B. Bramas, O. Coulaud
*
* This program runs the FMM Algorithm with the interpolation kernel based on
* uniform (grid points) interpolation (1/r kernel). It then compares the
* results with a direct computation.
*/
#include <string>
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
//
template<typename FReal, int ORDER>
using FInterpolationCell = FUnifCell<FReal, ORDER>;
template<typename FReal, typename GroupCellClass,
typename GroupContainerClass,
typename MatrixKernelClass, int ORDER>
using FInterpolationKernel = FUnifKernel<FReal,
GroupCellClass,
GroupContainerClass,
MatrixKernelClass,
ORDER> ;
const std::string interpolationType("Uniform interpolation");
#include "MPIInterpolationFMM.hpp"
......@@ -280,7 +280,7 @@ int main(int argc, char *argv[]) {
std::string outputFile(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "outputMPIFMA"));
FMpiFmaGenericWriter<FReal> paraWriter(outputFile,FMpiComm);
paraWriter.writeDistributionOfParticlesFromOctree(*computeOnGroupTree,loader.getNumberOfParticles(),mortonCellDistribution);
paraWriter.writeDistributionOfParticlesFromGroupedOctree(*computeOnGroupTree,loader.getNumberOfParticles(),mortonCellDistribution);
}
return 0;
}
......@@ -51,512 +51,520 @@
template <class FReal>
struct EwalParticle {
FPoint<FReal> position;
FReal forces[3];
FReal physicalValue;
FReal potential;
int index;
FPoint<FReal> position;
FReal forces[3];
FReal physicalValue;
FReal potential;
int index;
};
// Simply create particles and try the kernels
int main(int argc, char ** argv){
FHelpDescribeAndExit(argc, argv, "Please read the code to know more, sorry");
typedef double FReal;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
FHelpDescribeAndExit(argc, argv, "Please read the code to know more, sorry",
FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight,
FParameterDefinitions::InputFile,
FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads,
FParameterDefinitions::EnabledVerbose,
FParameterDefinitions::PeriodicityNbLevels);
typedef double FReal;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
#ifdef SCALFMM_USE_BLAS
// begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 13;
// typedefs
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebCell<FReal,ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
// begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 6;
// typedefs
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebCell<FReal,ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#else
typedef FSphericalCell<FReal> CellClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel<FReal, CellClass, ContainerClass > KernelClass;
const int DevP = FParameters::getValue(argc,argv,"-P", 9);
typedef FSphericalCell<FReal> CellClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel<FReal, CellClass, ContainerClass > KernelClass;
const int DevP = FParameters::getValue(argc,argv,"-P", 9);
#endif
using FmmClass = FFmmAlgorithmPeriodic<FReal,OctreeClass, CellClass, ContainerClass,
KernelClass, LeafClass > ;
using FmmClassNoPer= FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > ;
///////////////////////What we do/////////////////////////////
if( FParameters::existParameter(argc, argv, "-help")){
std::cout << ">> This executable has to be used to compute direct interaction either for periodic or non periodic system.\n";
std::cout << ">> options are -h H -sh SH [-per perdeep, -noper] -fin filenameIN (-bin) -[no]scale \n";
std::cout << ">> Recommended files : ../Data/EwalTest_Periodic.run ../Data/EwalTest_NoPeriodic.run\n";
std::cout << " Options " << std::endl;
std::cout << " -per perDeep " << std::endl;
std::cout << " -noper no periodic boundary conditions " << std::endl;
std::cout << " -verbose : print index x y z fx fy fy Q and V" << std::endl;
std::cout << " -noscale no scaling and we don't remove the dipole term " << std::endl;
exit(-1);
} //////////////////////////////////////////////////////////////
const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 4);
const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 2);
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 3);
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/EwalTest_Periodic.run");
// file for -saveError option
std::ofstream errorfile("outputEwaldError.txt",std::ios::out);
FTic counter;
// c conversion factor for coulombic terms in internal units
// c i.e. (unit(charge)**2/(4 pi eps0 unit(length))/unit(energy)
const FReal r4pie0 = FReal(138935.4835);
FReal scaleEnergy, scaleForce ;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
std::cout << "Opening : " << filename << "\n";
FDlpolyLoader<FReal> *loader = nullptr ;
if(FParameters::existParameter(argc, argv, "-bin")){
loader = new FDlpolyBinLoader<FReal>(filename);
}
else {
loader = new FDlpolyAsciiLoader<FReal>(filename);
}
if(! loader->isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
// ---------------------------------------------------
// DL_POLY CONSTANT
// ---------------------------------------------------
bool scale = true ;
if(FParameters::existParameter(argc, argv, "-noscale")){
scale = false ;
scaleEnergy = 1.0; // kcal mol^{-1}
scaleForce = 1.0 ; // 10 J mol^{-1} A^{-1}
}
else {
scaleEnergy = r4pie0 / 418.4 ; // kcal mol^{-1}
scaleForce = -r4pie0 ; // 10 J mol^{-1} A^{-1}
}
//
using FmmClass = FFmmAlgorithmPeriodic<FReal,OctreeClass, CellClass, ContainerClass,
KernelClass, LeafClass > ;
using FmmClassNoPer= FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > ;
///////////////////////What we do/////////////////////////////
if( FParameters::existParameter(argc, argv, "-help")){
std::cout << ">> This executable has to be used to compute direct interaction either for periodic or non periodic system.\n";
std::cout << ">> options are -h H -sh SH [-per perdeep, -noper] -fin filenameIN (-bin) -[no]scale \n";
std::cout << ">> Recommended files : ../Data/EwalTest_Periodic.run ../Data/EwalTest_NoPeriodic.run\n";
std::cout << " Options " << std::endl;
std::cout << " -per perDeep " << std::endl;
std::cout << " -noper no periodic boundary conditions " << std::endl;
std::cout << " -verbose : print index x y z fx fy fy Q and V" << std::endl;
std::cout << " -noscale no scaling and we don't remove the dipole term " << std::endl;
exit(-1);
} //////////////////////////////////////////////////////////////
const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 4);
const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 2);
const int PeriodicDeep = FParameters::getValue(argc,argv,FParameterDefinitions::PeriodicityNbLevels.options, 3);
const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/EwalTest_Periodic.run");
// file for -saveError option
std::ofstream errorfile("outputEwaldError.txt",std::ios::out);
FTic counter;
// c conversion factor for coulombic terms in internal units
// c i.e. (unit(charge)**2/(4 pi eps0 unit(length))/unit(energy)
const FReal r4pie0 = FReal(138935.4835);
FReal scaleEnergy, scaleForce ;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
std::cout << "Opening : " << filename << "\n";
FDlpolyLoader<FReal> *loader = nullptr ;
if( filename.find(".bin") != std::string::npos ){
loader = new FDlpolyBinLoader<FReal>(filename.c_str());
}
else {
loader = new FDlpolyAsciiLoader<FReal>(filename.c_str());
}
if(! loader->isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
// ---------------------------------------------------
// DL_POLY CONSTANT
// ---------------------------------------------------
bool scale = true ;
if(FParameters::existParameter(argc, argv, "-noscale")){
scale = false ;
scaleEnergy = 1.0; // kcal mol^{-1}
scaleForce = 1.0 ; // 10 J mol^{-1} A^{-1}
}
else {
scaleEnergy = r4pie0 / 418.4 ; // kcal mol^{-1}
scaleForce = -r4pie0 ; // 10 J mol^{-1} A^{-1}
}
//
#ifndef SCALFMM_USE_BLAS
CellClass::Init(DevP);
std::cout << " $$$$$$$$$$$$$$$ SPHERICAL VERSION $$$$$$$$$$$$"<<std::endl;
std::cout << " $$$$$$$$$$$$$$$ Order "<< DevP <<" $$$$$$$$$$$$"<<std::endl;
CellClass::Init(DevP);
std::cout << " $$$$$$$$$$$$$$$ SPHERICAL VERSION $$$$$$$$$$$$"<<std::endl;
std::cout << " $$$$$$$$$$$$$$$ Order "<< DevP <<" $$$$$$$$$$$$"<<std::endl;
#else
std::cout << " $$$$$$$$$$$$$$$ CHEBYCHEV VERSION $$$$$$$$$$$$" <<std::endl;
std::cout << " $$$$$$$$$$$$$$$ Order "<<ORDER <<" $$$$$$$$$$$$"<<std::endl;
std::cout << " $$$$$$$$$$$$$$$ CHEBYCHEV VERSION $$$$$$$$$$$$" <<std::endl;
std::cout << " $$$$$$$$$$$$$$$ Order "<<ORDER <<" $$$$$$$$$$$$"<<std::endl;
#endif
OctreeClass tree(NbLevels, SizeSubLevels, loader->getBoxWidth(), loader->getCenterOfBox());
// ---------------------------------------------------------------------------------
// Insert particles in the Octree
// --------------------------------------------------------------------------------- std::cout << "Creating & Inserting " << loader->getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX()
<< " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
FPoint<FReal> electricMoment(0.0,0.0,0.0) ;
EwalParticle<FReal> * const