Commit 862ef4b0 authored by berenger-bramas's avatar berenger-bramas

PreCalculs FMB

Threaded >> Thread

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@10 2616d619-271b-44dc-8df4-d4a8f33a7222
parent b2c72a49
......@@ -20,3 +20,4 @@ tmp/
*~
*TAGS
*#*#
......@@ -425,7 +425,7 @@ public:
* @return true if we can move down
*/
bool canProgressToDown() const {
return !this->isAtLeafLevel();
return !isAtLeafLevel();
}
/**
......@@ -433,7 +433,7 @@ public:
* @return true if we are at the bottom of the tree
*/
bool isAtLeafLevel() const {
return this->level() + 1 == OctreeHeight;
return level() + 1 == OctreeHeight;
}
/**
......@@ -491,17 +491,15 @@ public:
// To be able to access octree root
friend class Iterator;
///////////////////////////////////////////////////////////////////////////
// This part is related to the FMM algorithm (needed by M2M,M2L,etc.)
///////////////////////////////////////////////////////////////////////////
/** This function return a cell (if it exists) from a morton index and a level
* @param inIndex the index of the desired cell
* @param inLevel the level of the desired cell (cannot be infered from the index)
* @return the cell if it exist or null (0)
* This function starts from the root until it finds a missing cell or the right cell
* This function starts from the root until it find a missing cell or the right cell
*/
CellClass* getCell(const MortonIndex inIndex, const int inLevel){
SubOctreeTypes workingTree;
......@@ -570,8 +568,7 @@ public:
}
/** This function return an adresse of cell array from a morton index and a level
* So after calling this function you can test 8 cells (for each case if not null
* there is a cell allocated)
*
* @param inIndex the index of the desired cell array has to contains
* @param inLevel the level of the desired cell (cannot be infered from the index)
* @return the cell if it exist or null (0)
......@@ -603,19 +600,7 @@ public:
return &workingTree.tree->cellsAt(levelInTree)[treeLeafMask & inIndex];
}
/** This function fill a array with the distant neighbors of a cell
* Distant neighbors are child of the father neighbor (if not a direct
* neighbors)
* There is a maximum of 8 * 26 (3*3*3-1) distant neighbors
* @param inNeighbors the array to store the elements
* @param inIndex the index of the element we want the neighbors
* @param inLevel the level of the element
* @return the number of neighbors
*
* First find the direct neighbors of the cell (without brother)
* Second then for each neighbors of the parent of the cell
* Add child of this neighbors if they are not in the direct neighbors
*/
int getDistantNeighbors(CellClass* inNeighbors[208], const MortonIndex inIndex, const int inLevel){
// Take the neighbors != brothers
CellClass* directNeighbors[26];
......@@ -675,10 +660,10 @@ public:
}
/** This function return the particules list (leaf) from a morton index
/** This function return a cell (if it exists) from a morton index and a level
* @param inIndex the index of the desired cell
* @param inLevel the level of the desired cell (cannot be infered from the index)
* @return the list if it exists or null (0)
* @return the cell if it exist or null (0)
*
*/
FList<ParticuleClass*>* getLeaf(const MortonIndex inIndex){
......@@ -700,16 +685,17 @@ public:
return workingTree.leafTree->getLeaf(treeLeafMask & inIndex);
}
/** This function fill an array with the neighbors of a leaf
/** This function fill an array with the neighbors of a cell
* @param inNeighbors the array to store the elements
* @param inIndex the index of the element we want the neighbors
* @param inLevel the level of the element
* @return the number of neighbors
*/
int getLeafsNeighbors(FList<ParticuleClass*>* inNeighbors[26], const MortonIndex inIndex){
int getLeafsNeighbors(FList<ParticuleClass*>* inNeighbors[26], const MortonIndex inIndex, const int inLevel){
FTreeCoordinate center;
center.setPositionFromMorton(inIndex, leafIndex);
center.setPositionFromMorton(inIndex, inLevel);
const long limite = FMath::pow(2,leafIndex);
const long limite = FMath::pow(2,inLevel);
int idxNeighbors = 0;
......@@ -726,7 +712,7 @@ public:
// if we are not on the current cell
if( !(!idxX && !idxY && !idxZ) ){
const FTreeCoordinate other(center.getX() + idxX,center.getY() + idxY,center.getZ() + idxZ);
const MortonIndex mortonOther = other.getMortonIndex(leafIndex);
const MortonIndex mortonOther = other.getMortonIndex(inLevel);
// get cell
FList<ParticuleClass*>* const leaf = getLeaf(mortonOther);
// add to list if not null
......
......@@ -160,27 +160,27 @@ public:
/** Suboctree height accessor (leaf level + 1)
* @return subOctreeHeight */
int getSubOctreeHeight() const{
return this->subOctreeHeight;
return subOctreeHeight;
}
/** Suboctree position in the real tree
* @return subOctreePosition */
int getSubOctreePosition() const {
return this->subOctreePosition;
return subOctreePosition;
}
/** Return the more left leaf index
* the smallest index on the leafs array
* @return leftLeafIndex */
long getLeftLeafIndex() const {
return this->leftLeafIndex;
return leftLeafIndex;
}
/** Return the more right leaf index
* the biggest index on the leafs array
* @return rightLeafIndex */
long getRightLeafIndex() const {
return this->rightLeafIndex;
return rightLeafIndex;
}
/** Return the array of cells at a specious index
......@@ -188,19 +188,19 @@ public:
* @return cells[level] */
CellClass** cellsAt(const int level){
assert(level < subOctreeHeight, "Level out of memory", __LINE__, __FILE__);
return this->cells[level];
return cells[level];
}
/** To know if it is the root suboctree
* @return true if has parent otherwise return false */
bool hasParent() const {
return this->parent;
return parent;
}
/** To get access to the parent suboctree
* @return parent */
FAbstractSubOctree* getParent(){
return this->parent;
return parent;
}
/** To get the index of the current suboctree in the parent leafs array
......@@ -208,7 +208,7 @@ public:
* in this suboctree
* @return indexInParent */
long getIndexInParent() const{
return this->indexInParent;
return indexInParent;
}
};
......
......@@ -38,7 +38,7 @@ public:
* @return the position of the current cell
*/
MortonIndex getMortonIndex() const {
return this->index;
return index;
}
/**
......@@ -46,7 +46,7 @@ public:
* @param inPos the position given by the basic loader
*/
void setMortonIndex(const MortonIndex inIndex) {
this->index = inIndex;
index = inIndex;
}
};
......
......@@ -46,7 +46,7 @@ public:
* @return the position of the current cell
*/
F3DPosition getPosition() const {
return this->pos;
return pos;
}
/**
......@@ -54,7 +54,7 @@ public:
* @param inPos the position given by the basic loader
*/
void setPosition(const F3DPosition& inPos) {
this->pos = inPos;
pos = inPos;
}
};
......
......@@ -84,7 +84,7 @@ public:
// We need the current cell and the child
// child is an array (of 8 child) that may be null
kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild());
} while( octreeIterator.moveRight() );
} while(octreeIterator.moveRight());
octreeIterator.moveUp();
octreeIterator.gotoLeft();
}
......@@ -103,10 +103,10 @@ public:
// for each levels
for(int idxLevel = 2 ; idxLevel < OctreeHeight ; ++idxLevel ){
// for each cells
do {
const int numberOfNeighbors = tree->getDistantNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),idxLevel);
kernels->M2L( octreeIterator.getCurrentCell() , neighbors, numberOfNeighbors);
} while( octreeIterator.moveRight() );
do{
const int counter = tree->getDistantNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),idxLevel);
kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter);
} while(octreeIterator.moveRight());
octreeIterator.gotoLeft();
octreeIterator.moveDown();
}
......@@ -141,8 +141,8 @@ public:
do{
kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentList());
// need the current particules and neighbors particules
const int numberOfNeighbors = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex());
kernels->P2P( octreeIterator.getCurrentList() , neighbors, numberOfNeighbors);
const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),OctreeHeight-1);
kernels->P2P( octreeIterator.getCurrentList() , neighbors, counter);
} while(octreeIterator.moveRight());
FDEBUG( FDebug::Controller.write("Finished\n"); )
......
This diff is collapsed.
......@@ -64,7 +64,7 @@ public:
* Default destructor, simply close the file
*/
virtual ~FBasicLoader(){
this->file.close();
file.close();
}
/**
......
......@@ -26,8 +26,8 @@ private:
* Called by Fapplication instance to set the current app
* @param inCurrentApp current app
*/
static void SetCurrentApp(FAbstractApplication* const inCurrentApp){
CurrentApp = inCurrentApp;
static void SetCurrentApp(FAbstractApplication* const inCurrentApp){
CurrentApp = inCurrentApp;
}
/** To set CurrentApp */
......
#ifndef FCOMPLEXE_HPP
#define FCOMPLEXE_HPP
// /!\ Please, you must read the license at the bottom of this page
#include "FMath.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class
* Please read the license
*
* Propose basic complexe class
*/
class FComplexe {
double imag;
double real;
public:
FComplexe()
: imag(0), real(0){
}
FComplexe(const double inImag, const double inReal)
: imag(inImag), real(inReal){
}
FComplexe(const FComplexe& other)
: imag(other.imag), real(other.real){
}
FComplexe& operator=(const FComplexe& other){
this->imag = other.imag;
this->real = other.real;
return *this;
}
bool operator==(const FComplexe& other){
return FMath::LookEqual(this->imag,other.imag)
&& FMath::LookEqual(this->real,other.real);
}
bool operator!=(const FComplexe& other){
return !(*this == other);
}
double getImag() const{
return this->imag;
}
double getReal() const{
return this->real;
}
void setImag(const double inImag) {
this->imag = inImag;
}
void setReal(const double inReal) {
this->real = inReal;
}
};
#endif //FCOMPLEXE_HPP
// [--LICENSE--]
......@@ -32,7 +32,8 @@ private:
std::ostream* stream; //< Standart c++ ostream
/** Default constructor forbiden */
FDebug() : stream(&std::cout) {
FDebug(){
stream = &std::cout;
}
/** Default destructor forbiden */
......@@ -46,8 +47,8 @@ private:
* after this call stream is useless
*/
void close(){
this->stream->flush();
if(this->stream != &std::cout) delete(this->stream);
stream->flush();
if(stream != &std::cout) delete(stream);
}
/**
......@@ -77,7 +78,7 @@ public:
std::ofstream* const file = new std::ofstream();
file->open(filename);
this->stream = file;
stream = file;
}
/**
......@@ -85,7 +86,7 @@ public:
*/
void writeToCout(){
close();
this->stream = &std::cout;
stream = &std::cout;
}
/**
......@@ -95,7 +96,7 @@ public:
*/
template <class T>
FDebug& operator<<(const T& inMessage){
(*this->stream) << inMessage;
(*stream) << inMessage;
return *this;
}
......@@ -106,7 +107,7 @@ public:
*/
template <class T>
FDebug& write(const T& inMessage){
(*this->stream) << inMessage;
(*stream) << inMessage;
return *this;
}
......@@ -127,7 +128,7 @@ public:
oss << "Message from " << inFilePosition << " (at line " << inLinePosition <<")\n";
oss << ">> " << inMessage << "\n";
(*this->stream) << oss.str();
(*stream) << oss.str();
return *this;
}
......@@ -148,7 +149,7 @@ public:
std::ostringstream oss;
oss << "[Value] " << inVariable << " = " << inValue << " at line " << inLinePosition <<" (file " << inFilePosition << ")\n";
(*this->stream) << oss.str();
(*stream) << oss.str();
return *this;
}
......
......@@ -12,6 +12,9 @@
* Propose basic math functions or indirections
*/
struct FMath{
static const double FPi;
static const double FPiDiv2;
/** To get absolute value */
template <class NumType>
static NumType Abs(const NumType inV){
......@@ -51,8 +54,22 @@ struct FMath{
static bool Between(const NumType inValue, const NumType inMin, const NumType inMax){
return ( inMin <= inValue && inValue < inMax ? true : false);
}
/** To get sqrt of a double */
static double Sqrt(const double inValue){
return sqrt(inValue);
}
/** To get sqrt of a double */
static double Atan2(const double inValue1,const double inValue2){
return atan2(inValue1,inValue2);
}
/** To get sqrt of a double */
static double Sin(const double inValue){
return sin(inValue);
}
};
const double FMath::FPi = 3.14159265358979323846;
const double FMath::FPiDiv2 = 1.57079632679489661923;
#endif //FMATH_HPP
......
......@@ -5,6 +5,8 @@
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include "../Sources/Containers/FOctree.hpp"
#include "../Sources/Containers/FList.hpp"
......@@ -17,10 +19,8 @@
#include "../Sources/Core/FFMMAlgorithm.hpp"
#include "../Sources/Core/FSimpleKernels.hpp"
#include "../Sources/Utils/FTic.hpp"
// We use openmp to count time (portable and easy to manage)
// Compile by : g++ testFMMAlgorithm.cpp ../Sources/Utils/FDebug.cpp ../Sources/Utils/FAssertable.cpp -O2 -o testFMMAlgorithm.exe
// Compile by : g++ testFMMAlgorithm.cpp ../Sources/Utils/FAssertable.cpp -lgomp -fopenmp -O2 -o testFMMAlgorithm.exe
/** This program show an example of use of
* the fmm basic algo
......@@ -145,44 +145,43 @@ public:
int main(int , char ** ){
const int NbLevels = 10;//10;
const int NbSubLevels = 3;//3
const long NbPart = 200000;//2E6;
const long NbPart = 20000;//2E6;
MyTestParticule* particules = new MyTestParticule[NbPart];
FTic counter;
srand ( 1 ); // volontary set seed to constant
// -----------------------------------------------------
std::cout << "Creating " << NbPart << " particules ..." << std::endl;
counter.tic();
const double CreatingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
new (&particules[idxPart]) MyTestParticule(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
const double CreatingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (CreatingEndTime-CreatingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
FOctree<MyTestParticule, MyTestCell, NbLevels, NbSubLevels> tree(1.0,F3DPosition(0.5,0.5,0.5));
// -----------------------------------------------------
std::cout << "Inserting particules ..." << std::endl;
counter.tic();
const double InsertingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
tree.insert(&particules[idxPart]);
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
const double InsertingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (InsertingEndTime-InsertingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Working on particules ..." << std::endl;
counter.tic();
const double WorkingStartTime = omp_get_wtime();
MyTestKernels<MyTestParticule, MyTestCell> kernels;
FFMMAlgorithm<MyTestParticule, MyTestCell, NbLevels, NbSubLevels> algo(&tree,&kernels);
algo.execute();
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
const double WorkingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (WorkingEndTime-WorkingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
......@@ -217,13 +216,13 @@ int main(int , char ** ){
// -----------------------------------------------------
std::cout << "Deleting particules ..." << std::endl;
counter.tic();
const double DeletingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
particules[idxPart].~MyTestParticule();
}
delete [] particules;
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
const double DeletingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (DeletingEndTime-DeletingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
return 0;
......
// /!\ Please, you must read the license at the bottom of this page
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include "../Sources/Containers/FOctree.hpp"
#include "../Sources/Containers/FList.hpp"
#include "../Sources/Utils/FAssertable.hpp"
#include "../Sources/Utils/F3DPosition.hpp"
#include "../Sources/Core/FBasicParticule.hpp"
#include "../Sources/Core/FBasicCell.hpp"
#include "../Sources/Core/FFMMAlgorithm.hpp"
#include "../Sources/Core/FFmbKernels.hpp"
// We use openmp to count time (portable and easy to manage)
// Compile by : g++ testFMMAlgorithm.cpp ../Sources/Utils/FAssertable.cpp -lgomp -fopenmp -O2 -o testFMMAlgorithm.exe
/** This program show an example of use of
* the fmm basic algo
* it also check that eachh particules is little or longer
* related that each other
*/
/** Custom particule class */
class MyTestParticule : public FBasicParticule{
protected:
// To store data during downard pass
long dataDown;
public:
MyTestParticule(const double inX, const double inY, const double inZ) : FBasicParticule(inX,inY,inZ) {
this->dataDown = 0;
}
MyTestParticule(const F3DPosition& inPos) : FBasicParticule(inPos) {
this->dataDown = 0;
}
MyTestParticule(){
this->dataDown = 0;
}
long getDataDown() const {
return this->dataDown;
}
void setDataDown(const long inData){
this->dataDown = inData;
}
};
/** Custom cell */
class MyTestCell : public FBasicCell {
// To store data during upward and downward pass
long dataUp, dataDown;
public:
MyTestCell(const MortonIndex inIndex) : FBasicCell(inIndex) {
this->dataUp = 0;
this->dataDown = 0;
}
MyTestCell(){
this->dataUp = 0;
this->dataDown = 0;
}
long getDataUp() const {
return this->dataUp;
}
void setDataUp(const long inData){
this->dataUp = inData;
}
long getDataDown() const {
return this->dataDown;
}
void setDataDown(const long inData){
this->dataDown = inData;
}
};
/** Custom Kernel
* This kernel simply store in each element the number of elements
* it represents
*/
template< class ParticuleClass, class CellClass>
class MyTestKernels : public FAbstractKernels<ParticuleClass,CellClass> {
public:
// Before upward
virtual void P2M(CellClass* const pole, FList<ParticuleClass*>* const particules) {
// the pole represents all particules under
pole->setDataUp(particules->getSize());
}
// During upward
virtual void M2M(CellClass* const pole, CellClass** const child) {
// A parent represents the sum of the child
for(int idx = 0 ; idx < 8 ; ++idx){
if(child[idx]){
pole->setDataUp(pole->getDataUp() + child[idx]->getDataUp());
}
}
}
// Before Downward
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size) {
// The pole is impacted by what represent other poles
for(int idx = 0 ; idx < size ; ++idx){
pole->setDataDown(pole->getDataDown() + distantNeighbors[idx]->getDataUp());
}
}
// During Downward
virtual void L2L(CellClass* const local, CellClass** const child) {
// Each child is impacted by the father
for(int idx = 0 ; idx < 8 ; ++idx){
if(child[idx]){
child[idx]->setDataDown(local->getDataDown() + child[idx]->getDataDown());
}
}
}
// After Downward
virtual void L2P(CellClass* const local, FList<ParticuleClass*>* const particules){
// The particules is impacted by the parent cell
typename FList<ParticuleClass*>::BasicIterator iter(*particules);
while( iter.isValide() ){
iter.value()->setDataDown(iter.value()->getDataDown() + local->getDataDown());
iter.progress();
}
}
// After Downward
virtual void P2P(FList<ParticuleClass*>* const currentBox, FList<ParticuleClass*>** directNeighbors, const int size) {
// Each particules targeted is impacted by the particules sources
long inc = currentBox->getSize() - 1;
for(int idx = 0 ; idx < size ; ++idx){
inc += directNeighbors[idx]->getSize();
}
typename FList<ParticuleClass*>::BasicIterator iter(*currentBox);
while( iter.isValide() ){
iter.value()->setDataDown(iter.value()->getDataDown() + inc);
iter.progress();
}
}
};
// Simply create particules and try the kernels
int main(int , char ** ){
const int NbLevels = 10;//10;
const int NbSubLevels = 3;//3
const long NbPart = 20000;//2E6;
const int BoxWidth = 2;//1
MyTestParticule* particules = new MyTestParticule[NbPart];
srand ( 1 ); // volontary set seed to constant
// -----------------------------------------------------
std::cout << "Creating " << NbPart << " particules ..." << std::endl;
const double CreatingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
new (&particules[idxPart]) MyTestParticule(double(rand())/RAND_MAX,double(rand())/RAND_MAX,double(rand())/RAND_MAX);
}
const double CreatingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (CreatingEndTime-CreatingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
FOctree<MyTestParticule, MyTestCell, NbLevels, NbSubLevels> tree(BoxWidth,F3DPosition(0.0,0.0,0.0));
// -----------------------------------------------------
std::cout << "Inserting particules ..." << std::endl;
const double InsertingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
tree.insert(&particules[idxPart]);
}
const double InsertingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (InsertingEndTime-InsertingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Working on particules ..." << std::endl;
const double WorkingStartTime = omp_get_wtime();
FFmbKernels<MyTestParticule, MyTestCell> kernels(NbLevels,BoxWidth);
//FFMMAlgorithm<MyTestParticule, MyTestCell, NbLevels, NbSubLevels> algo(&tree,&kernels);
//algo.execute();
const double WorkingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (WorkingEndTime-WorkingStartTime) << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Deleting particules ..." << std::endl;
const double DeletingStartTime = omp_get_wtime();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
particules[idxPart].~MyTestParticule();
}
delete [] particules;
const double DeletingEndTime = omp_get_wtime();
std::cout << "Done " << "(" << (DeletingEndTime-DeletingStartTime) << "s)." << std::endl;
// -----------------------------------------------------