Commit 7ad5f4c5 authored by berenger-bramas's avatar berenger-bramas
Browse files

Fmb Kernels

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@11 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 862ef4b0
......@@ -106,7 +106,7 @@ public:
*/
void insert(ParticuleClass* const inParticule){
const MortonIndex particuleIndex = getLeafMortonFromPosition(inParticule->getPosition());
root.insert( particuleIndex, inParticule, height);
root.insert( particuleIndex, inParticule, height, this->boxWidthAtLevel);
}
......
......@@ -5,6 +5,7 @@
#include "../Utils/F3DPosition.hpp"
#include "../Utils/FAssertable.hpp"
#include "../Utils/FMath.hpp"
#include "../Utils/FConvert.hpp"
#include "FTreeCoordinate.hpp"
#include "FList.hpp"
......@@ -71,12 +72,17 @@ protected:
* when computing
* @param arrayIndex the index at the leaf index of the new element
*/
void createPreviousCells(MortonIndex arrayIndex, MortonIndex inLeafCellIndex){
void createPreviousCells(MortonIndex arrayIndex, MortonIndex inLeafCellIndex, const double* const inBoxWidthAtLevel){
int indexLevel = this->subOctreeHeight - 1;
while(indexLevel >= 0 && !this->cells[indexLevel][arrayIndex]){
this->cells[indexLevel][arrayIndex] = new CellClass();
this->cells[indexLevel][arrayIndex]->setMortonIndex(inLeafCellIndex);
const int realLevel = indexLevel + this->getSubOctreePosition();
this->cells[indexLevel][arrayIndex]->setPosition(FConvert::MortonToPosition(inLeafCellIndex,realLevel,inBoxWidthAtLevel[realLevel]));
--indexLevel;
arrayIndex >>= 3;
inLeafCellIndex >>= 3;
}
......@@ -87,8 +93,8 @@ protected:
* for example it updates the leaf array marges and calls createPreviousCells()
* @param arrayIndex the position of the new leaf in the leafs array
*/
void newLeafInserted(const long arrayIndex, const MortonIndex inLeafCellIndex){
createPreviousCells(arrayIndex,inLeafCellIndex);
void newLeafInserted(const long arrayIndex, const MortonIndex inLeafCellIndex, const double* const inBoxWidthAtLevel){
createPreviousCells(arrayIndex,inLeafCellIndex, inBoxWidthAtLevel);
// Update if this is the bottom left
if(arrayIndex < this->leftLeafIndex) this->leftLeafIndex = arrayIndex;
if(arrayIndex > this->rightLeafIndex) this->rightLeafIndex = arrayIndex;
......@@ -151,7 +157,7 @@ public:
* @param inParticule the particule to insert (must inherite from FAbstractParticule)
* @param inParticule the inTreeHeight the height of the tree
*/
virtual void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight) = 0;
virtual void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight, const double* const inBoxWidthAtLevel) = 0;
///////////////////////////////////////
// This is the FOctree::Iterator Part
......@@ -273,13 +279,13 @@ public:
/**
* Refer to FAbstractSubOctree::insert
*/
void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight){
void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight, const double* const inBoxWidthAtLevel){
// Get the morton index for the leaf level
const MortonIndex arrayIndex = FAbstractSubOctree<ParticuleClass,CellClass>::getLeafIndex(index,inTreeHeight);
// is there already a leaf?
if( !this->leafs[arrayIndex] ){
this->leafs[arrayIndex] = new FList<ParticuleClass*>();
FAbstractSubOctree<ParticuleClass,CellClass>::newLeafInserted( arrayIndex , index);
FAbstractSubOctree<ParticuleClass,CellClass>::newLeafInserted( arrayIndex , index, inBoxWidthAtLevel);
}
// add particule to leaf list
this->leafs[arrayIndex]->pushFront(inParticule);
......@@ -353,7 +359,7 @@ public:
/**
* Refer to FAbstractSubOctree::insert
*/
void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight){
void insert(const MortonIndex index, ParticuleClass* const inParticule, const int inTreeHeight, const double* const inBoxWidthAtLevel){
// We need the morton index at the bottom level of this sub octree
// so we remove the right side
const MortonIndex arrayIndex = FAbstractSubOctree<ParticuleClass,CellClass>::getLeafIndex(index,inTreeHeight);
......@@ -373,10 +379,10 @@ public:
}
// We need to inform parent class
FAbstractSubOctree<ParticuleClass,CellClass>::newLeafInserted( arrayIndex, index >> (3 * (inTreeHeight-nextSubOctreePosition) ) );
FAbstractSubOctree<ParticuleClass,CellClass>::newLeafInserted( arrayIndex, index >> (3 * (inTreeHeight-nextSubOctreePosition) ), inBoxWidthAtLevel);
}
// Ask next suboctree to insert the particule
this->subleafs[arrayIndex]->insert( index, inParticule, inTreeHeight);
this->subleafs[arrayIndex]->insert( index, inParticule, inTreeHeight, inBoxWidthAtLevel);
}
/** To get access to leafs elements (child suboctree)
......
......@@ -31,6 +31,12 @@ public:
* @param inIndex the position of the current cell
*/
virtual void setMortonIndex(const MortonIndex inIndex) = 0;
/**
* Must be implemented by each user Cell class
* @param inPosition the position of the current cell
*/
virtual void setPosition(const F3DPosition& inPosition) = 0;
};
......
......@@ -23,13 +23,13 @@ public:
virtual void P2M(CellClass* const pole, FList<ParticuleClass*>* const particules) = 0;
/** M2M */
virtual void M2M(CellClass* const pole, CellClass** const child) = 0;
virtual void M2M(CellClass* const pole, CellClass** const child, const int inLevel) = 0;
/** M2L */
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size) = 0;
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size, const int inLevel) = 0;
/** L2L */
virtual void L2L(CellClass* const pole, CellClass** const child) = 0;
virtual void L2L(CellClass* const pole, CellClass** const child, const int inLevel) = 0;
/** L2P */
virtual void L2P(CellClass* const pole, FList<ParticuleClass*>* const particules) = 0;
......
......@@ -15,18 +15,11 @@
class FBasicCell : public FAbstractCell {
protected:
MortonIndex index; //< Cell's position in the tree
F3DPosition position; //< Cell's spacial position
public:
/**
* Constructor with a position
* @param inIndex the morton index
*/
FBasicCell(const MortonIndex inIndex) : index(inIndex) {
}
/** Default constructor */
FBasicCell(){
FBasicCell() : index(0), position(0,0,0){
}
/** Default destructor */
......@@ -38,7 +31,7 @@ public:
* @return the position of the current cell
*/
MortonIndex getMortonIndex() const {
return index;
return this->index;
}
/**
......@@ -46,7 +39,21 @@ public:
* @param inPos the position given by the basic loader
*/
void setMortonIndex(const MortonIndex inIndex) {
index = inIndex;
this->index = inIndex;
}
/**
* @return the position of the current cell
*/
F3DPosition getPosition() const{
return this->position;
}
/**
* @param inPosition the position of the current cell
*/
void setPosition(const F3DPosition& inPosition){
this->position = inPosition;
}
};
......
......@@ -4,11 +4,13 @@
#include "../Utils/FAssertable.hpp"
#include "../Utils/FDebug.hpp"
#include "../Utils/FTic.hpp"
#include "../Containers/FOctree.hpp"
#include "FAbstractKernels.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFMMAlgorithm
......@@ -23,6 +25,8 @@ class FFMMAlgorithm : public FAssertable{
FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>* const tree; //< The octree to work on
FAbstractKernels<ParticuleClass, CellClass>* const kernels; //< The kernels
FDEBUG(FTic counter);
public:
/** The constructor need the octree and the kernels used for computation
* @param inTree the octree
......@@ -56,7 +60,8 @@ public:
/** P2M */
void bottomPass(){
FDEBUG( FDebug::Controller.write("Start Bottom Pass\n") );
FDEBUG( FDebug::Controller.write("\tStart Bottom Pass\n").write(FDebug::Flush) );
FDEBUG(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
// Iterate on leafs
......@@ -67,12 +72,14 @@ public:
kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentList());
} while(octreeIterator.moveRight());
FDEBUG( FDebug::Controller.write("Finished\n"); )
FDEBUG(counter.tac(););
FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
}
/** M2M */
void upwardPass(){
FDEBUG( FDebug::Controller.write("Start Upward Pass\n"); );
FDEBUG( FDebug::Controller.write("\tStart Upward Pass\n").write(FDebug::Flush); );
FDEBUG(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
......@@ -83,18 +90,20 @@ public:
do{
// 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());
kernels->M2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), idxLevel);
} while(octreeIterator.moveRight());
octreeIterator.moveUp();
octreeIterator.gotoLeft();
}
FDEBUG( FDebug::Controller.write("Finished\n"); )
FDEBUG(counter.tac(););
FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
}
/** M2L L2L */
void downardPass(){
FDEBUG( FDebug::Controller.write("Start Downward Pass\n"); );
FDEBUG( FDebug::Controller.write("\tStart Downward Pass (M2L)\n").write(FDebug::Flush); );
FDEBUG(counter.tic(););
{ // first M2L
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
......@@ -105,12 +114,17 @@ public:
// for each cells
do{
const int counter = tree->getDistantNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),idxLevel);
kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter);
kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, idxLevel);
} while(octreeIterator.moveRight());
octreeIterator.gotoLeft();
octreeIterator.moveDown();
}
}
FDEBUG(counter.tac(););
FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
FDEBUG( FDebug::Controller.write("\tStart Downward Pass (L2L)\n").write(FDebug::Flush); );
FDEBUG(counter.tic(););
{ // second L2L
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.moveDown();
......@@ -119,19 +133,22 @@ public:
for(int idxLevel = 2 ; idxLevel < heightMinusOne ; ++idxLevel ){
// for each cells
do{
kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild());
kernels->L2L( octreeIterator.getCurrentCell() , octreeIterator.getCurrentChild(), idxLevel);
} while(octreeIterator.moveRight());
octreeIterator.gotoLeft();
octreeIterator.moveDown();
}
}
FDEBUG( FDebug::Controller.write("Finished\n"); )
FDEBUG(counter.tac(););
FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
}
/** P2P */
void directPass(){
FDEBUG( FDebug::Controller.write("Start Direct Pass\n"); );
FDEBUG( FDebug::Controller.write("\tStart Direct Pass\n").write(FDebug::Flush); );
FDEBUG(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
......@@ -145,7 +162,8 @@ public:
kernels->P2P( octreeIterator.getCurrentList() , neighbors, counter);
} while(octreeIterator.moveRight());
FDEBUG( FDebug::Controller.write("Finished\n"); )
FDEBUG(counter.tac(););
FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
}
......
This diff is collapsed.
......@@ -31,17 +31,17 @@ public:
}
/** Print the morton index */
virtual void M2M(CellClass* const pole, CellClass** const child) {
virtual void M2M(CellClass* const pole, CellClass** const child, const int inLevel) {
FDEBUG( FDebug::Controller << "M2M : " << pole->getMortonIndex() << "\n" );
}
/** Print the morton index */
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size) {
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size, const int inLevel) {
FDEBUG( FDebug::Controller << "M2L : " << pole->getMortonIndex() << " (" << size << ")\n" );
}
/** Print the morton index */
virtual void L2L(CellClass* const pole, CellClass** const child) {
virtual void L2L(CellClass* const pole, CellClass** const child, const int inLevel) {
FDEBUG( FDebug::Controller << "L2L : " << pole->getMortonIndex() << "\n" );
}
......
......@@ -5,6 +5,7 @@
// To get memcpy
#include <cstring>
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class F3DPosition
......@@ -131,6 +132,30 @@ public:
return *this;
}
/**
* Subtract to all dim the other position
* @param other the value to substract
* @return the current object after being subtracted
*/
F3DPosition& operator-=(const F3DPosition& other){
this->x -= other.x;
this->y -= other.y;
this->z -= other.z;
return *this;
}
/**
* Affect to all dim the other position
* @param other the value to afect
* @return the current object after being affected
*/
F3DPosition& operator+=(const F3DPosition& other){
this->x += other.x;
this->y += other.y;
this->z += other.z;
return *this;
}
};
/**
......@@ -159,6 +184,32 @@ F3DPosition operator+(const F3DPosition& inPosition, const double inValue){
return position;
}
/**
* Operator F3Position minus F3Position
* This substract one from anther
* @param inPosition the position to reduce
* @param inOther the position to decrease/substract inPosition
* @return the resulting position
*/
F3DPosition operator-(const F3DPosition& inPosition, const F3DPosition& inOther){
F3DPosition position(inPosition);
position -= inOther;
return position;
}
/**
* Operator F3Position plus F3Position
* This substract one from anther
* @param inPosition the position to reduce
* @param inOther the position to increase inPosition
* @return the resulting position
*/
F3DPosition operator+(const F3DPosition& inPosition, const F3DPosition& inOther){
F3DPosition position(inPosition);
position += inOther;
return position;
}
#endif //F3DPOSITION_HPP
// [--LICENSE--]
......@@ -52,6 +52,30 @@ public:
void setReal(const double inReal) {
this->real = inReal;
}
FComplexe& operator+=(const FComplexe& other){
this->real += other.real;
this->imag += other.imag;
return *this;
}
void mulRealAndImag(const double inValue){
this->imag *= inValue;
this->real *= inValue;
}
void inc(const double inIncReal, const double inIncImag){
this->real += inIncReal;
this->imag += inIncImag;
}
void incReal(const double inIncReal){
this->real += inIncReal;
}
void incImag(const double inIncImag){
this->imag += inIncImag;
}
};
......
#ifndef FCONVERT_HPP
#define FCONVERT_HPP
#include "../Containers/FTreeCoordinate.hpp"
#include "../Utils/F3DPosition.hpp"
/** This class proposes some convert functions
*
*/
class FConvert {
public :
/** To get spatial position (F3DPosition) from morton data
* @param inIndex the morton index of the object
* @param inLevel the level of the current object
* @param inWidthAtLevel the width of the box at this level
* return outPosition the result
*/
static F3DPosition MortonToPosition(const MortonIndex inIndex, const int inLevel, const double inWidthAtLevel){
FTreeCoordinate treePosition;
treePosition.setPositionFromMorton(inIndex, inLevel);
F3DPosition outPosition(
treePosition.getX() * inWidthAtLevel + inWidthAtLevel/2,
treePosition.getY() * inWidthAtLevel + inWidthAtLevel/2,
treePosition.getZ() * inWidthAtLevel + inWidthAtLevel/2
);
return outPosition;
}
};
#endif
......@@ -33,7 +33,7 @@ private:
/** Default constructor forbiden */
FDebug(){
stream = &std::cout;
this->stream = &std::cout;
}
/** Default destructor forbiden */
......@@ -47,8 +47,8 @@ private:
* after this call stream is useless
*/
void close(){
stream->flush();
if(stream != &std::cout) delete(stream);
flush();
if(this->stream != &std::cout) delete(this->stream);
}
/**
......@@ -78,7 +78,7 @@ public:
std::ofstream* const file = new std::ofstream();
file->open(filename);
stream = file;
this->stream = file;
}
/**
......@@ -86,7 +86,7 @@ public:
*/
void writeToCout(){
close();
stream = &std::cout;
this->stream = &std::cout;
}
/**
......@@ -96,8 +96,7 @@ public:
*/
template <class T>
FDebug& operator<<(const T& inMessage){
(*stream) << inMessage;
return *this;
return write(inMessage);
}
/**
......@@ -107,10 +106,31 @@ public:
*/
template <class T>
FDebug& write(const T& inMessage){
(*stream) << inMessage;
(*this->stream) << inMessage;
return *this;
}
/** Flush data into stream */
void flush(){
this->stream->flush();
}
enum FlushType{
Flush,
FlushWithLine,
};
/**
* stream operator to flush debug data
* @param inType flush type
* @return current FDebug
*/
FDebug& write(const FlushType inType){
if(inType == FlushWithLine) (*this->stream) << '\n';
flush();
return *this;
}
/**
* to write debug data with line & file
* @param inMessage a message - from any type - to print
......@@ -128,7 +148,7 @@ public:
oss << "Message from " << inFilePosition << " (at line " << inLinePosition <<")\n";
oss << ">> " << inMessage << "\n";
(*stream) << oss.str();
(*this->stream) << oss.str();
return *this;
}
......@@ -149,7 +169,7 @@ public:
std::ostringstream oss;
oss << "[Value] " << inVariable << " = " << inValue << " at line " << inLinePosition <<" (file " << inFilePosition << ")\n";
(*stream) << oss.str();
(*this->stream) << oss.str();
return *this;
}
......
......@@ -19,6 +19,9 @@ int main(void){
// Write a developer information
FDEBUG( FDebug::Controller.writeFromLine("Strange things are there!", __LINE__, __FILE__); )
// Flush
FDEBUG( FDebug::Controller << FDebug::Flush );
// Change stream type
FDEBUG( FDebug::Controller.writeToFile("testDebug.out.temp"); )
FDEBUG( FDebug::Controller << "Hello Wordl 2 the return\n");
......
......@@ -5,7 +5,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include "../Sources/Utils/FTic.hpp"
#include "../Sources/Containers/FOctree.hpp"
#include "../Sources/Containers/FList.hpp"
......@@ -19,8 +19,7 @@
#include "../Sources/Core/FFMMAlgorithm.hpp"
#include "../Sources/Core/FSimpleKernels.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
// Compile by : g++ testFMMAlgorithm.cpp ../Sources/Utils/FAssertable.cpp ../Sources/Utils/FDebug.cpp -O2 -o testFMMAlgorithm.exe
/** This program show an example of use of
* the fmm basic algo
......@@ -57,10 +56,6 @@ 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;
......@@ -92,7 +87,7 @@ public:
pole->setDataUp(particules->getSize());
}
// During upward
virtual void M2M(CellClass* const pole, CellClass** const child) {
virtual void M2M(CellClass* const pole, CellClass** const child, const int inLevel) {
// A parent represents the sum of the child
for(int idx = 0 ; idx < 8 ; ++idx){
if(child[idx]){
......@@ -101,14 +96,14 @@ public:
}
}
// Before Downward
virtual void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size) {