diff --git a/Sources/Core/FAbstractCell.hpp b/Sources/Core/FAbstractCell.hpp
index 8b418bfefffb9094c19de775690284466e549ecf..80cc533509ae9945d7da73bb3767e33247113d26 100644
--- a/Sources/Core/FAbstractCell.hpp
+++ b/Sources/Core/FAbstractCell.hpp
@@ -11,6 +11,12 @@
*
* This class define the method that every cell class
* has to implement.
+*
+* In fact FOctree & FFMMAlgorithm need this function to be implemented.
+* But you cannot use this interface with the extension (as an example :
+* because the compiler will faill to know if getMortonIndex is coming
+* from this interface or from the extension)
+*
* @warning Inherite from this class when implement a specific cell type
*/
class FAbstractCell{
diff --git a/Sources/Core/FAbstractParticule.hpp b/Sources/Core/FAbstractParticule.hpp
index 5b0034253fabf30c3b935bdafd0b7734f6eaa13d..803fe390c68fffa772c45afb217f0876dfa905b7 100644
--- a/Sources/Core/FAbstractParticule.hpp
+++ b/Sources/Core/FAbstractParticule.hpp
@@ -13,6 +13,13 @@ class F3DPosition;
*
* This class define the method that every particule class
* has to implement.
+*
+* In fact FOctree & FFMMAlgorithm need this function to be implemented.
+* But you cannot use this interface with the extension (as an example :
+* because the compiler will faill to know if getPosition is coming
+* from this interface or from the extension)
+*
+*
* @warning Inherite from this class when implement a specific particule type
*/
class FAbstractParticule{
diff --git a/Sources/Core/FBasicCell.hpp b/Sources/Core/FBasicCell.hpp
index 2ab33ad6ae33c1ccef0c32004abe5f4794e760cd..09d8401f1a978707c7bcb5a0dd868970b166b432 100644
--- a/Sources/Core/FBasicCell.hpp
+++ b/Sources/Core/FBasicCell.hpp
@@ -2,59 +2,25 @@
#define FBASICCELL_HPP
// /!\ Please, you must read the license at the bottom of this page
-#include "FAbstractCell.hpp"
-#include "../Containers/FTreeCoordinate.hpp"
+#include "FExtendPosition.hpp"
+#include "FExtendMortonIndex.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FBasicCell
* Please read the license
*
-* This class defines a basic cell used for examples.
+* This class defines a basic cell used for examples. It extends
+* the mininum, only what is needed by FOctree and FFMMAlgorithm
+* to make the things working.
+* By using this extension it will implement the FAbstractCell without
+* inheriting from it.
*/
-class FBasicCell : public FAbstractCell {
-protected:
- MortonIndex index; //< Cell's position in the tree
- F3DPosition position; //< Cell's spacial position
-
+class FBasicCell : public FExtendPosition, public FExtendMortonIndex {
public:
- /** Default constructor */
- FBasicCell() : index(0), position(0,0,0){
- }
-
/** Default destructor */
virtual ~FBasicCell(){
}
-
- /**
- * From the FAbstractParticule definition
- * @return the position of the current cell
- */
- MortonIndex getMortonIndex() const {
- return this->index;
- }
-
- /**
- * This function is needed by the basic loader to fill the current particule
- * @param inPos the position given by the basic loader
- */
- void setMortonIndex(const MortonIndex 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;
- }
};
diff --git a/Sources/Core/FBasicKernels.hpp b/Sources/Core/FBasicKernels.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2601bf256714403ab33e02a488965e4757e23844
--- /dev/null
+++ b/Sources/Core/FBasicKernels.hpp
@@ -0,0 +1,62 @@
+#ifndef FBASICKERNELS_HPP
+#define FBASICKERNELS_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include "FAbstractKernels.hpp"
+
+#include "../Utils/FDebug.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class AbstractKernels
+* @brief
+* Please read the license
+*
+* This kernels simply shows the details of the information
+* it receives (in debug)
+*/
+template< class ParticuleClass, class CellClass>
+class FBasicKernels : public FAbstractKernels<ParticuleClass,CellClass> {
+public:
+ /** Default destructor */
+ virtual ~FBasicKernels(){
+ }
+
+ /** When init the kernel */
+ virtual void init(){}
+
+ /** Print the number of particules */
+ virtual void P2M(CellClass* const pole, FList<ParticuleClass*>* const particules) {
+ FDEBUG( FDebug::Controller << "P2M : " << particules->getSize() << "\n" );
+ }
+
+ /** Print the morton index */
+ 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, const int inLevel) {
+ FDEBUG( FDebug::Controller << "M2L : " << pole->getMortonIndex() << " (" << size << ")\n" );
+ }
+
+ /** Print the morton index */
+ virtual void L2L(CellClass* const pole, CellClass** const child, const int inLevel) {
+ FDEBUG( FDebug::Controller << "L2L : " << pole->getMortonIndex() << "\n" );
+ }
+
+ /** Print the number of particules */
+ virtual void L2P(CellClass* const pole, FList<ParticuleClass*>* const particules){
+ FDEBUG( FDebug::Controller << "L2P : " << particules->getSize() << "\n" );
+ }
+
+ /** Print the number of particules */
+ virtual void P2P(FList<ParticuleClass*>* const currentBox, FList<ParticuleClass*>** directNeighbors, const int size) {
+ FDEBUG( FDebug::Controller << "P2P : " << currentBox->getSize() << " (" << size << ")\n" );
+ }
+};
+
+
+#endif //FBASICKERNELS_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FBasicParticule.hpp b/Sources/Core/FBasicParticule.hpp
index c6306939712627479c951ef7ede0c5f36d94cbb6..d6692d7a145b7750c29a30676a0f89b4bfec98d0 100644
--- a/Sources/Core/FBasicParticule.hpp
+++ b/Sources/Core/FBasicParticule.hpp
@@ -2,60 +2,24 @@
#define FBASICPARTICULE_HPP
// /!\ Please, you must read the license at the bottom of this page
-#include "FAbstractParticule.hpp"
-#include "../Utils/F3DPosition.hpp"
+#include "FExtendPosition.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FBasicParticule
* Please read the license
*
-* This class defines a basic particule used for examples.
+* This class defines a basic particule used for examples. It extends
+* the mininum, only what is needed by FOctree and FFMMAlgorithm
+* to make the things working.
+* By using this extension it will implement the FAbstractParticule without
+* inheriting from it.
*/
-class FBasicParticule : public FAbstractParticule{
-protected:
- F3DPosition pos; //< Particule's position
-
+class FBasicParticule : public FExtendPosition{
public:
- /**
- * Constructor with a position
- * @param inX x position
- * @param inY y position
- * @param inZ z position
- */
- FBasicParticule(const double inX, const double inY, const double inZ) : pos(inX,inY,inZ) {
- }
-
- /**
- * Constructor with a position object
- * @param inPos particule position
- */
- FBasicParticule(const F3DPosition& inPos) : pos(inPos) {
- }
-
- /** Default constructor */
- FBasicParticule(){
- }
-
/** Default destructor */
virtual ~FBasicParticule(){
}
-
- /**
- * From the FAbstractParticule definition
- * @return the position of the current cell
- */
- F3DPosition getPosition() const {
- return pos;
- }
-
- /**
- * This function is needed by the basic loader to fill the current particule
- * @param inPos the position given by the basic loader
- */
- void setPosition(const F3DPosition& inPos) {
- pos = inPos;
- }
};
diff --git a/Sources/Core/FExtendForces.hpp b/Sources/Core/FExtendForces.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bacdfbe1772403c9bf0d303a770859fddf8015ac
--- /dev/null
+++ b/Sources/Core/FExtendForces.hpp
@@ -0,0 +1,60 @@
+#ifndef FEXTENDFORCES_HPP
+#define FEXTENDFORCES_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include "../Utils/F3DPosition.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendForces
+* Please read the license
+*
+* This class is an extenssion.
+* It proposes a 3d array as a forces vector.
+*/
+class FExtendForces {
+protected:
+ F3DPosition forces; //< 3D vector stored in a position object
+
+public:
+ /** Default constructor */
+ FExtendForces() {
+ }
+
+ /** Copy constructor */
+ FExtendForces(const FExtendForces& other) : forces(other.forces) {
+ }
+
+ /** Destructor */
+ virtual ~FExtendForces(){
+ }
+
+ /** Copy operator */
+ FExtendForces& operator=(const FExtendForces& other) {
+ this->forces = other.forces;
+ return *this;
+ }
+
+ /** Return the forces */
+ F3DPosition getForces() const {
+ return this->forces;
+ }
+
+ /** Set Forces */
+ void setForces(const F3DPosition& inForces) {
+ this->forces = inForces;
+ }
+
+ /** Set Forces with 3 doubles */
+ void setForces(const double inFx, const double inFy, const double inFz) {
+ this->forces.setX(inFx);
+ this->forces.setY(inFy);
+ this->forces.setZ(inFz);
+ }
+
+};
+
+
+#endif //FEXTENDFORCES_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FExtendMortonIndex.hpp b/Sources/Core/FExtendMortonIndex.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7e7d5c0523f01a4de1d819bf72051e029a1ad315
--- /dev/null
+++ b/Sources/Core/FExtendMortonIndex.hpp
@@ -0,0 +1,52 @@
+#ifndef FEXTENDMORTONINDEX_HPP
+#define FEXTENDMORTONINDEX_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include "../Containers/FTreeCoordinate.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendMortonIndex
+* Please read the license
+* This class is an extenssion.
+* It proposes a mortonIndex.
+*/
+class FExtendMortonIndex {
+protected:
+ MortonIndex mortonIndex; //< Morton index (need by most elements)
+
+public:
+ /** Default constructor */
+ FExtendMortonIndex() : mortonIndex(0) {
+ }
+
+ /** Copy constructor */
+ FExtendMortonIndex(const FExtendMortonIndex& other) : mortonIndex(other.mortonIndex) {
+ }
+
+ /** Destructor */
+ virtual ~FExtendMortonIndex(){
+ }
+
+ /** Copy operator */
+ FExtendMortonIndex& operator=(const FExtendMortonIndex& other) {
+ this->mortonIndex = other.mortonIndex;
+ return *this;
+ }
+
+ /** To get the morton index */
+ MortonIndex getMortonIndex() const {
+ return this->mortonIndex;
+ }
+
+ /** To set the morton index */
+ void setMortonIndex(const MortonIndex inMortonIndex) {
+ this->mortonIndex = inMortonIndex;
+ }
+
+};
+
+
+#endif //FEXTENDMORTONINDEX_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FExtendPosition.hpp b/Sources/Core/FExtendPosition.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3688dec3fd3e48700bb971877ebc95b1d8fb1d44
--- /dev/null
+++ b/Sources/Core/FExtendPosition.hpp
@@ -0,0 +1,59 @@
+#ifndef FEXTENDPOSITION_HPP
+#define FEXTENDPOSITION_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include "../Utils/F3DPosition.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendPosition
+* Please read the license
+* This class is an extenssion.
+* It proposes a mortonIndex.
+*/
+class FExtendPosition {
+protected:
+ F3DPosition position; //< The position
+
+public:
+ /** Default constructor */
+ FExtendPosition() {
+ }
+
+ /** Copy constructor */
+ FExtendPosition(const FExtendPosition& other) : position(other.position) {
+ }
+
+ /** Destructor */
+ virtual ~FExtendPosition(){
+ }
+
+ /** Copy operator */
+ FExtendPosition& operator=(const FExtendPosition& other) {
+ this->position = other.position;
+ return *this;
+ }
+
+ /** To get the position */
+ F3DPosition getPosition() const {
+ return this->position;
+ }
+
+ /** To set the position */
+ void setPosition(const F3DPosition& inPosition) {
+ this->position = inPosition;
+ }
+
+ /** To set the position from 3 doubles */
+ void setPosition(const double inX, const double inY, const double inZ) {
+ this->position.setX(inX);
+ this->position.setY(inY);
+ this->position.setZ(inZ);
+ }
+
+};
+
+
+#endif //FEXTENDPOSITION_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FExtendPotential.hpp b/Sources/Core/FExtendPotential.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..65f272af9a8cd88b1401c8ce120720c75d8fe896
--- /dev/null
+++ b/Sources/Core/FExtendPotential.hpp
@@ -0,0 +1,50 @@
+#ifndef FEXTENDPOTENTIAL_HPP
+#define FEXTENDPOTENTIAL_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendPotential
+* Please read the license
+* This class is an extenssion.
+* It proposes a Potential (double).
+*/
+class FExtendPotential {
+protected:
+ double potential; //< The potential extended
+
+public:
+ /** Default constructor */
+ FExtendPotential() : potential(0) {
+ }
+
+ /** Copy constructor */
+ FExtendPotential(const FExtendPotential& other) : potential(other.potential) {
+ }
+
+ /** Destructor */
+ virtual ~FExtendPotential(){
+ }
+
+ /** Copy operator */
+ FExtendPotential& operator=(const FExtendPotential& other) {
+ this->potential = other.potential;
+ return *this;
+ }
+
+ /** To get the potential */
+ double getPotential() const {
+ return this->potential;
+ }
+
+ /** To set the potential */
+ void setPotential(const double inPotential) {
+ this->potential = inPotential;
+ }
+
+};
+
+
+#endif //FEXTENDPOTENTIAL_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FExtendValue.hpp b/Sources/Core/FExtendValue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f6aaf8ca81eea31218d32061322997da25d62145
--- /dev/null
+++ b/Sources/Core/FExtendValue.hpp
@@ -0,0 +1,50 @@
+#ifndef FEXTENDVALUE_HPP
+#define FEXTENDVALUE_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendValue
+* Please read the license
+* This class is an extenssion.
+* It proposes a value (double).
+*/
+class FExtendValue {
+protected:
+ double value; //< A simple value
+
+public:
+ /** Default constructor */
+ FExtendValue() : value(0) {
+ }
+
+ /** Copy constructor */
+ FExtendValue(const FExtendValue& other) : value(other.value) {
+ }
+
+ /** Destructor */
+ virtual ~FExtendValue(){
+ }
+
+ /** Copy Constructor */
+ FExtendValue& operator=(const FExtendValue& other) {
+ this->value = other.value;
+ return *this;
+ }
+
+ /** To get the value */
+ double getValue() const {
+ return this->value;
+ }
+
+ /** To set the value */
+ void setValue(const double invalue) {
+ this->value = invalue;
+ }
+
+};
+
+
+#endif //FEXTENDVALUE_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Core/FFMMAlgorithm.hpp b/Sources/Core/FFMMAlgorithm.hpp
index d0ede7273e625e78d5d90938609c414187a6c196..cf712c6972469d7f4f31e77324b0ccf6878d4d9a 100644
--- a/Sources/Core/FFMMAlgorithm.hpp
+++ b/Sources/Core/FFMMAlgorithm.hpp
@@ -25,7 +25,7 @@ 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);
+ FDEBUG_TIME(FTic counter); //< In case of debug count the time
public:
/** The constructor need the octree and the kernels used for computation
@@ -60,8 +60,8 @@ public:
/** P2M */
void bottomPass(){
- FDEBUG( FDebug::Controller.write("\tStart Bottom Pass\n").write(FDebug::Flush) );
- FDEBUG(counter.tic(););
+ FDEBUG_TRACE( FDebug::Controller.write("\tStart Bottom Pass\n").write(FDebug::Flush) );
+ FDEBUG_TIME(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
// Iterate on leafs
@@ -72,14 +72,14 @@ public:
kernels->P2M( octreeIterator.getCurrentCell() , octreeIterator.getCurrentList());
} while(octreeIterator.moveRight());
- FDEBUG(counter.tac(););
- FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
+ FDEBUG_TIME(counter.tac(););
+ FDEBUG_TRACE( FDebug::Controller << "\tFinished (") FDEBUG_TIME(<< counter.elapsed() <<) FDEBUG_TRACE("s)\n"; )
}
/** M2M */
void upwardPass(){
- FDEBUG( FDebug::Controller.write("\tStart Upward Pass\n").write(FDebug::Flush); );
- FDEBUG(counter.tic(););
+ FDEBUG_TRACE( FDebug::Controller.write("\tStart Upward Pass\n").write(FDebug::Flush); );
+ FDEBUG_TIME(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
@@ -96,14 +96,14 @@ public:
octreeIterator.gotoLeft();
}
- FDEBUG(counter.tac(););
- FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
+ FDEBUG_TIME(counter.tac(););
+ FDEBUG_TRACE( FDebug::Controller << "\tFinished (") FDEBUG_TIME(<< counter.elapsed() <<) FDEBUG_TRACE("s)\n"; )
}
/** M2L L2L */
void downardPass(){
- FDEBUG( FDebug::Controller.write("\tStart Downward Pass (M2L)\n").write(FDebug::Flush); );
- FDEBUG(counter.tic(););
+ FDEBUG_TRACE( FDebug::Controller.write("\tStart Downward Pass (M2L)\n").write(FDebug::Flush); );
+ FDEBUG_TIME(counter.tic(););
{ // first M2L
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
@@ -120,11 +120,11 @@ public:
octreeIterator.moveDown();
}
}
- FDEBUG(counter.tac(););
- FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
+ FDEBUG_TIME(counter.tac(););
+ FDEBUG_TRACE( FDebug::Controller << "\tFinished (") FDEBUG_TIME(<< counter.elapsed() <<) FDEBUG_TRACE("s)\n"; )
- FDEBUG( FDebug::Controller.write("\tStart Downward Pass (L2L)\n").write(FDebug::Flush); );
- FDEBUG(counter.tic(););
+ FDEBUG_TRACE( FDebug::Controller.write("\tStart Downward Pass (L2L)\n").write(FDebug::Flush); );
+ FDEBUG_TIME(counter.tic(););
{ // second L2L
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.moveDown();
@@ -140,15 +140,15 @@ public:
}
}
- FDEBUG(counter.tac(););
- FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
+ FDEBUG_TIME(counter.tac(););
+ FDEBUG_TRACE( FDebug::Controller << "\tFinished (") FDEBUG_TIME(<< counter.elapsed() <<) FDEBUG_TRACE("s)\n"; )
}
/** P2P */
void directPass(){
- FDEBUG( FDebug::Controller.write("\tStart Direct Pass\n").write(FDebug::Flush); );
- FDEBUG(counter.tic(););
+ FDEBUG_TRACE( FDebug::Controller.write("\tStart Direct Pass\n").write(FDebug::Flush); );
+ FDEBUG_TIME(counter.tic(););
typename FOctree<ParticuleClass, CellClass, OctreeHeight, SubtreeHeight>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
@@ -162,8 +162,8 @@ public:
kernels->P2P( octreeIterator.getCurrentList() , neighbors, counter);
} while(octreeIterator.moveRight());
- FDEBUG(counter.tac(););
- FDEBUG( FDebug::Controller << "\tFinished (" << counter.elapsed() << "s)\n"; )
+ FDEBUG_TIME(counter.tac(););
+ FDEBUG_TRACE( FDebug::Controller << "\tFinished (") FDEBUG_TIME(<< counter.elapsed() <<) FDEBUG_TRACE("s)\n"; )
}
diff --git a/Sources/Core/FFmaParticule.hpp b/Sources/Core/FFmaParticule.hpp
index d78e93acfa7f31d1f4ad43c0a8bcf8fcbb3aa67f..41d47a4ae9333092fdbe3c0e04aad813cb215531 100644
--- a/Sources/Core/FFmaParticule.hpp
+++ b/Sources/Core/FFmaParticule.hpp
@@ -3,60 +3,21 @@
// /!\ Please, you must read the license at the bottom of this page
#include "FBasicParticule.hpp"
+#include "FExtendValue.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFmaParticule
* Please read the license
*
-* This class defines a basic particule used for examples.
+* This class defines a particule for FMA loader.
+* As defined in FFmaLoader it needs {FBasicParticule,FExtendValue}
*/
-class FFmaParticule : public FBasicParticule{
-protected:
- double value; //< Particule's fma data
-
+class FFmaParticule : public FBasicParticule, public FExtendValue {
public:
- /**
- * Constructor with a position
- * @param inX x position
- * @param inY y position
- * @param inZ z position
- * @param inValue particule value
- */
- FFmaParticule(const double inX, const double inY, const double inZ, const double inValue)
- : FBasicParticule(inX,inY,inZ), value(inValue) {
- }
-
- /**
- * Constructor with a position object
- * @param inPos particule position
- */
- FFmaParticule(const F3DPosition& inPos) : FBasicParticule(inPos), value(0){
- }
-
- /** Default constructor */
- FFmaParticule() : value(0) {
- }
-
/** Default destructor */
virtual ~FFmaParticule(){
}
-
- /**
- * From the Fmb needed
- * @return the value of the current cell
- */
- double getValue() const {
- return this->value;
- }
-
- /**
- * This function is needed by the fma loader to fill the current particule
- * @param inValue the data given by the fma loader
- */
- void setValue(const double inValue) {
- this->value = inValue;
- }
};
diff --git a/Sources/Fmb/FExtendFmbCell.hpp b/Sources/Fmb/FExtendFmbCell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c682ca5768a8188a99f0c2f5db3e9410326fa98
--- /dev/null
+++ b/Sources/Fmb/FExtendFmbCell.hpp
@@ -0,0 +1,58 @@
+#ifndef FEXTENDFMBCELL_HPP
+#define FEXTENDFMBCELL_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include <string.h>
+
+#include "../Utils/FComplexe.hpp"
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class FExtendFmbCell
+* Please read the license
+*/
+template <int P>
+class FExtendFmbCell {
+protected:
+ static const int FMB_Info_P = P; //< P >> FMB_Info.P
+ static const int MultipoleSize = int(((FMB_Info_P)+1) * ((FMB_Info_P)+2) * 0.5); //< The size of the multipole
+
+ FComplexe multipole_exp[MultipoleSize]; //< For multipole extenssion
+ FComplexe local_exp[MultipoleSize]; //< For local extenssion
+
+public:
+ /** Default constructor */
+ FExtendFmbCell() {
+ }
+
+ /** Constructor */
+ FExtendFmbCell(const FExtendFmbCell& other){
+ (*this) = other;
+ }
+
+ /** Default destructor */
+ virtual ~FExtendFmbCell(){
+ }
+
+ /** Copy constructor */
+ FExtendFmbCell& operator=(const FExtendFmbCell& other) {
+ memcpy(multipole_exp, other.multipole_exp, sizeof(multipole_exp));
+ memcpy(local_exp, other.local_exp, sizeof(local_exp));
+ return *this;
+ }
+
+ /** Get Multipole */
+ FComplexe* getMultipole() {
+ return this->multipole_exp;
+ }
+ /** Get Local */
+ FComplexe* getLocal() {
+ return this->local_exp;
+ }
+
+};
+
+
+#endif //FEXTENDFMBCELL_HPP
+
+// [--LICENSE--]
diff --git a/Sources/Fmb/FFmbKernels.hpp b/Sources/Fmb/FFmbKernels.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce91004a350adea4640c460d7338cbaddb78be5e
--- /dev/null
+++ b/Sources/Fmb/FFmbKernels.hpp
@@ -0,0 +1,1160 @@
+#ifndef FFMBKERNELS_HPP
+#define FFMBKERNELS_HPP
+// /!\ Please, you must read the license at the bottom of this page
+
+#include "../Core/FAbstractKernels.hpp"
+
+#include "../Containers/FTreeCoordinate.hpp"
+#include "../Containers/FList.hpp"
+
+#include "../Utils/F3DPosition.hpp"
+#include "../Utils/FComplexe.hpp"
+#include "../Utils/FMath.hpp"
+
+#include <iostream>
+
+
+/**
+* @author Berenger Bramas (berenger.bramas@inria.fr)
+* @class AbstractKernels
+* @brief
+* Please read the license
+*
+* This kernels simply shows the details of the information
+* it receives
+*/
+template< class ParticuleClass, class CellClass>
+class FFmbKernels : public FAbstractKernels<ParticuleClass,CellClass> {
+ // P is a input parameter
+ static const int FMB_Info_P = 2;
+
+ //#ifdef _GRAVITATIONAL_
+ static const int FMB_Info_eps_soft_square = 1;
+ // else = 0
+ //#endif
+
+ // LMax is used in several algorithm
+ // it is just a copy of FMB_Info.P
+ static const int LMax = FMB_Info_P;
+
+ // Can be 2 * FMB_Info_P if user ask to
+ static const int FMB_Info_M2L_P = FMB_Info_P;
+ static const int FMB_Info_M2L_exp_size = ((FMB_Info_M2L_P)+1) * ((FMB_Info_M2L_P)+2) * 0.5;
+
+ // Default value set in main
+ static const int FMB_Info_ws = 1;
+
+ // INTERACTION_LIST_SIZE_ALONG_1_DIM
+ static const int size1Dim = (2*(2*(FMB_Info_ws)+1) +1);
+ // HALF_INTERACTION_LIST_SIZE_ALONG_1_DIM
+ static const int halphSize1Dim = (2*(FMB_Info_ws)+1);
+
+ // EXPANSION_SIZE(FMB_Info.P)
+ static const int FMB_Info_exp_size = ((FMB_Info_P)+1) * ((FMB_Info_P)+2) * 0.5;
+ // NEXP_SIZE(FMB_Info.P)
+ static const int FMB_Info_nexp_size = (FMB_Info_P + 1) * (FMB_Info_P + 1);
+
+ // Level of the tree
+ const int treeHeight;
+ // Width of the box at the root level
+ const double treeWidthAtRoot;
+
+ // transfer_M2M_container
+ FComplexe*** transitionM2M;
+ // transfer_L2L_container
+ FComplexe*** transitionL2L;
+
+ // transfer_container
+ FComplexe***** transferM2L;
+
+ //[OK] spherical_harmonic_Outer_coefficients_array
+ double sphereHarmoOuterCoef[FMB_Info_exp_size];
+ //[OK] spherical_harmonic_Inner_coefficients_array
+ double sphereHarmoInnerCoef[FMB_Info_exp_size];
+
+ FComplexe current_thread_Y[FMB_Info_exp_size];
+
+ // p_Y_theta_derivated
+ FComplexe current_thread_Y_theta_derivated[FMB_Info_exp_size];
+
+ // [OK?] p_precomputed_cos_and_sin_array
+ FComplexe cosSin[FMB_Info_M2L_P + 1];
+ // [OK?] p_associated_Legendre_function_Array
+ double legendre[FMB_Info_M2L_exp_size];
+
+ // pow_of_I_array
+ static const double PiArrayInner[4];
+ // pow_of_O_array
+ static const double PiArrayOuter[4];
+
+ // To store spherical position
+ struct Spherical {
+ double r, cosTheta, sinTheta, phi;
+ };
+
+ int expansion_Redirection_array_for_j[FMB_Info_M2L_P + 1 ];
+
+ F3DPosition force_sum;
+ double potential_sum;
+
+ //////////////////////////////////////////////////////////////////
+ // Allocation
+ //////////////////////////////////////////////////////////////////
+
+ void expansion_Redirection_array_for_j_Initialize() {
+ for( int h = 0; h <= FMB_Info_M2L_P ; ++h ){
+ expansion_Redirection_array_for_j[h] = static_cast<int>( h * ( h + 1 ) * 0.5 );
+ }
+ }
+
+ //spherical_harmonic_Outer_and_Inner_coefficients_array_Initialize
+ void sphericalHarmonicInitialize(){
+ // Outer coefficients:
+ //std::cout << "sphereHarmoOuterCoef\n";
+ double factOuter = 1.0;
+ // in FMB code stoped at <= FMB_Info_M2L_P but this is not sufficient
+ for(int idxP = 0 ; idxP < FMB_Info_exp_size; factOuter *= (++idxP) ){
+ this->sphereHarmoOuterCoef[idxP] = factOuter;
+ //std::cout << this->sphereHarmoOuterCoef[idxl] << "\n";
+ }
+
+ // Inner coefficients:
+ double* currentInner = this->sphereHarmoInnerCoef;
+ double factInner = 1.0;
+ double powN1idxP = 1.0;
+ //std::cout << "sphereHarmoInnerCoef\n";
+ for(int idxP = 0 ; idxP <= this->FMB_Info_M2L_P ; factInner *= (++idxP), powN1idxP = -powN1idxP){
+ for(int idxMP = 0, fact_l_m = factInner; idxMP <= idxP ; fact_l_m *= idxP+(++idxMP), ++currentInner){
+ *currentInner = powN1idxP / fact_l_m;
+ //std::cout << (*currentInner) << "\n";
+ }
+ }
+
+ //for(int temp = 0 ; temp < 6 ; ++temp){
+ // std::cout << this->sphereHarmoInnerCoef[temp] << "\n";
+ //}
+ }
+
+
+ // transfer_L2L_Allocate
+ // transfer_M2M_Allocate
+ // transfer_M2L_Allocate
+ void transferAllocate(){
+ // M2M L2L
+ this->transitionM2M = new FComplexe**[this->treeHeight];
+ this->transitionL2L = new FComplexe**[this->treeHeight];
+
+ for(int idxLevel = 0 ; idxLevel < this->treeHeight ; ++idxLevel){
+
+ this->transitionM2M[idxLevel] = new FComplexe*[8];
+ this->transitionL2L[idxLevel] = new FComplexe*[8];
+
+ for(long idxChild = 0; idxChild < 8; ++idxChild){
+ this->transitionM2M[idxLevel][idxChild] = new FComplexe[FMB_Info_exp_size];
+ this->transitionL2L[idxLevel][idxChild] = new FComplexe[FMB_Info_exp_size];
+ }
+ }
+ // M2L
+ this->transferM2L = new FComplexe****[this->treeHeight+1];
+
+ for(int idxLevel = 0; idxLevel < this->treeHeight; ++idxLevel){
+ this->transferM2L[idxLevel] = new FComplexe***[this->size1Dim];
+
+ for(long idxD1 = 0 ; idxD1 < this->size1Dim; ++idxD1){
+ this->transferM2L[idxLevel][idxD1] = new FComplexe**[this->size1Dim];
+
+ for(long idxD2 = 0; idxD2 < this->size1Dim; ++idxD2){
+ this->transferM2L[idxLevel][idxD1][idxD2] = new FComplexe*[this->size1Dim];
+
+ for(long idxD3 = 0; idxD3 < this->size1Dim; ++idxD3){
+ const int x = idxD1 - this->halphSize1Dim;
+ const int y = idxD2 - this->halphSize1Dim;
+ const int z = idxD3 - this->halphSize1Dim;
+
+ if( ( x*x + y*y + z*z ) >= ( 3*this->FMB_Info_ws*this->FMB_Info_ws + 0.1 ) ){
+ //[Blas] this->transferM2L[idxLevel][idxD1][idxD2][idxD3] = new FComplexe[this->FMB_Info_exp_size * this->FMB_Info_nexp_size];
+ this->transferM2L[idxLevel][idxD1][idxD2][idxD3] = new FComplexe[this->FMB_Info_M2L_exp_size];
+ }
+ else {
+ this->transferM2L[idxLevel][idxD1][idxD2][idxD3] = NULL;
+ }
+ }
+ }
+ }
+ }
+ }
+ // transfer_L2L_free
+ // transfer_M2M_free
+ // transfer_M2L_free
+ void transferDeallocate(){
+ // M2M L2L
+ for(long idxLevel = 0 ; idxLevel < this->treeHeight ; ++idxLevel){
+ for(long idxChild = 0; idxChild < 8; ++idxChild){
+ delete [] this->transitionM2M[idxLevel][idxChild];
+ delete [] this->transitionL2L[idxLevel][idxChild];
+ }
+ delete [] this->transitionM2M[idxLevel];
+ delete [] transitionL2L[idxLevel];
+ }
+ delete [] this->transitionM2M;
+ delete [] this->transitionL2L;
+ // M2L
+ for(int idxLevel = 0 ; idxLevel < this->treeHeight; ++idxLevel){
+ for(long idxD1 = 0 ; idxD1 < this->size1Dim ; ++idxD1){
+ for(long idxD2 = 0 ; idxD2 < this->size1Dim ; ++idxD2){
+ for(long idxD3 = 0 ; idxD3 < this->size1Dim; ++idxD3){
+ delete [] this->transferM2L[idxLevel][idxD1][idxD2][idxD3];
+ }
+ delete [] this->transferM2L[idxLevel][idxD1][idxD2];
+ }
+ delete [] this->transferM2L[idxLevel][idxD1];
+ }
+ delete [] this->transferM2L[idxLevel];
+ }
+ delete [] this->transferM2L;
+ }
+
+ //////////////////////////////////////////////////////////////////
+ // Utils
+ //////////////////////////////////////////////////////////////////
+
+ // position_2_r_cos_th_sin_th_ph
+ void positionToSphere(const F3DPosition& inVector, Spherical* const outSphere ){
+ const double x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY());
+
+ outSphere->r = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ()));
+ outSphere->phi = FMath::Atan2(inVector.getY(),inVector.getX());
+ outSphere->cosTheta = inVector.getZ() / outSphere->r; // cos_th = z/r
+ outSphere->sinTheta = FMath::Sqrt(x2y2) / outSphere->r; // sin_th = sqrt(x^2 + y^2)/r
+ }
+
+ // associated_Legendre_function_Fill_complete_array_of_values_for_cos
+ void legendreFunction( const double inCosTheta, const double inSinTheta, double* const outResults ){
+ // l=0: results[current++] = 1.0; // P_0^0(cosTheta) = 1
+ int idxCurrent = 0;
+ outResults[idxCurrent++] = 1.0;
+
+ // l=1:
+ // Compute P_1^{0} using (3) and store it into results_array:
+ outResults[idxCurrent++] = inCosTheta;
+
+ // Compute P_1^1 using (2 bis) and store it into results_array
+ const double invSinTheta = -inSinTheta; // somx2 => -sinTheta
+ outResults[idxCurrent++] = invSinTheta;
+
+ // l>1:
+ int idxCurrent1m = 1; //pointer on P_{l-1}^m P_1^0
+ int idxCurrent2m = 0; //pointer on P_{l-2}^m P_0^0
+ double fact = 3.0;
+
+ // Remark: p_results_array_l_minus_1_m and p_results_array_l_minus_2_m
+ // just need to be incremented at each iteration.
+ for(int idxl = 2; idxl <= this->LMax ; ++idxl ){
+ for( int idxm = 0; idxm <= idxl - 2 ; ++idxm , ++idxCurrent , ++idxCurrent1m , ++idxCurrent2m ){
+ // Compute P_l^m, l >= m+2, using (1) and store it into results_array:
+ outResults[idxCurrent] = (inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m] - ( idxl + idxm - 1 )
+ * outResults[idxCurrent2m] ) / ( idxl - idxm );
+ }
+ // p_results_array_l_minus_1_m now points on P_{l-1}^{l-1}
+
+ // Compute P_l^{l-1} using (3) and store it into ptrResults:
+ outResults[idxCurrent++] = inCosTheta * ( 2 * idxl - 1 ) * outResults[idxCurrent1m];
+
+ // Compute P_l^l using (2 bis) and store it into results_array:
+ outResults[idxCurrent++] = fact * invSinTheta * outResults[idxCurrent1m];
+
+ fact += 2.0;
+ ++idxCurrent1m;
+ }
+ }
+
+ // spherical_harmonic_Inner
+ void harmonicInner(const Spherical& inSphere, FComplexe* const outResults){
+
+ FComplexe* ptrCosSin = this->cosSin;
+ for(int idxl = 0 , idxlMod4 = 0; idxl <= LMax ; ++idxl, ++idxlMod4, ++ptrCosSin){
+ if(idxlMod4 == 4) idxlMod4 = 0;
+ const double angle = idxl * inSphere.phi + this->PiArrayInner[idxlMod4];
+
+ ptrCosSin->setReal( FMath::Sin(angle + FMath::FPiDiv2) );
+ ptrCosSin->setImag( FMath::Sin(angle) );
+ //std::cout<< idxl << "=" << ptrCosSin->getReal() << "/" << ptrCosSin->getImag() << " (" << idxl << "/" << inSphere.phi << "/" << PiArray[lMod4] << ")\n";
+ }
+
+ legendreFunction(inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+ //printf("FMB_Info_M2L_exp_size=%d\n",FMB_Info_M2L_exp_size);
+ //for(int temp = 0 ; temp < FMB_Info_M2L_exp_size ; ++temp){
+ // printf("%f\n",this->legendre[temp]);
+ //}
+
+ FComplexe* currentResult = outResults;
+ int idxLegendre = 0;//ptr_associated_Legendre_function_Array
+ int idxSphereHarmoCoef = 0;
+ double idxRl = 1.0 ;
+
+ for(int idxl = 0; idxl <= this->LMax ; ++idxl, idxRl *= inSphere.r){
+ for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxSphereHarmoCoef, ++idxLegendre){
+ const double magnitude = this->sphereHarmoInnerCoef[idxSphereHarmoCoef] * idxRl * legendre[idxLegendre];
+ currentResult->setReal( magnitude * this->cosSin[idxm].getReal() );
+ currentResult->setImag( magnitude * this->cosSin[idxm].getImag() );
+ //printf("magnitude=%f idxRl=%f sphereHarmoInnerCoef=%f real=%f imag=%f\n",magnitude,idxRl,this->sphereHarmoInnerCoef[idxSphereHarmoCoef],currentResult->getReal(),currentResult->getImag());
+ }
+ }
+
+ }
+ // spherical_harmonic_Outer
+ void harmonicOuter(const Spherical& inSphere, FComplexe* const outResults){
+
+ FComplexe* ptrCosSin = this->cosSin;
+ for(int idxl = 0, idxlMod4 = 0; idxl <= LMax ; ++idxl, ++idxlMod4, ++ptrCosSin){
+ if(idxlMod4 == 4) idxlMod4 = 0;
+ const double angle = idxl * inSphere.phi + this->PiArrayOuter[idxlMod4];
+
+ ptrCosSin->setReal( FMath::Sin(angle + FMath::FPiDiv2) );
+ ptrCosSin->setImag( FMath::Sin(angle) );
+
+ //printf("l=%d \t inSphere.phi=%f \t this->PiArrayOuter[idxlMod4]=%f \t angle=%f \t FMath::Sin(angle + FMath::FPiDiv2)=%f \t FMath::Sin(angle)=%f\n",
+ // idxl, inSphere.phi, this->PiArrayOuter[idxlMod4], angle, FMath::Sin(angle + FMath::FPiDiv2) , FMath::Sin(angle));
+ }
+
+ legendreFunction(inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+
+ int idxLegendre = 0;
+ FComplexe* currentResult = outResults;
+
+ const double invR = 1/inSphere.r;
+ double idxRl1 = invR;
+ for(int idxl = 0 ; idxl <= this->LMax ; ++idxl, idxRl1 *= invR){
+ for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxLegendre){
+ const double magnitude = this->sphereHarmoOuterCoef[idxl-idxm] * idxRl1 * this->legendre[idxLegendre];
+ currentResult->setReal( magnitude * this->cosSin[idxm].getReal() );
+ currentResult->setImag( magnitude * this->cosSin[idxm].getImag() );
+ //printf("l=%d\t m=%d\t idxRl1=%f\t magnitude=%f\n",idxl,idxm,idxRl1,magnitude);
+ //printf("l=%d\t m=%d\t this->cosSin[idxm].getReal()=%f\t this->cosSin[idxm].getImag()=%f\n",
+ // idxl,idxm,this->cosSin[idxm].getReal(),this->cosSin[idxm].getImag());
+ }
+ }
+ }
+
+ /** spherical_harmonic_Inner_and_theta_derivated
+ * Returns the value of the partial derivative of the spherical harmonic
+ *relative to theta. We have for all m such that -(l-1) <= m <= l-1 :
+ *
+ *(d H_l^m(theta, phi))/(d theta)
+ *= (-1)^m sqrt((l-|m|)!/(l+|m|)!) (d P_l^{|m|}(cos theta))/(d theta) e^{i.m.phi}
+ *= (-1)^m sqrt((l-|m|)!/(l+|m|)!) 1/sqrt{1-cos^2 theta} [l cos theta P_l^{|m|}(cos theta) - (l+|m|) P_{l-1}^{|m|}(cos theta) ] e^{i.m.phi}
+ *Since theta is in the range [0, Pi], we have: sin theta > 0 and therefore
+ *sqrt{1-cos^2 theta} = sin theta. Thus:
+ *
+ *(d H_l^m(theta, phi))/(d theta)
+ *= (-1)^m sqrt((l-|m|)!/(l+|m|)!) 1/(sin theta) [l cos theta P_l^{|m|}(cos theta) - (l+|m|) P_{l-1}^{|m|}(cos theta) ] e^{i.m.phi}
+ *For |m|=l, we have~:
+ *(d H_l^l(theta, phi))/(d theta)
+ *= (-1)^m sqrt(1/(2l)!) 1/(sin theta) [l cos theta P_l^l(cos theta) ] e^{i.m.phi}
+ *
+ *Remark: for 0<m<=l:
+ *(d H_l^{-m}(theta, phi))/(d theta) = [(d H_l^{-m}(theta, phi))/(d theta)]*
+ *
+ *
+ *
+ *Therefore, we have for (d Inner_l^m(r, theta, phi))/(d theta):
+ *
+ *|m|<l: (d Inner_l^m(r, theta, phi))/(d theta) =
+ *(i^m (-1)^l / (l+|m|)!) 1/(sin theta) [l cos theta P_l^{|m|}(cos theta) - (l+|m|) P_{l-1}^{|m|}(cos theta) ] e^{i.m.phi} r^l
+ *|m|=l: (d Inner_l^m(r, theta, phi))/(d theta) =
+ *(i^m (-1)^l / (l+|m|)!) 1/(sin theta) [l cos theta P_l^l(cos theta) ] e^{i.m.phi} r^l
+ *
+ *
+ */
+ void harmonicInnerThetaDerivated(
+ const Spherical& inSphere,
+ FComplexe * results_array,
+ FComplexe * theta_derivated_results_array
+ ){
+ FComplexe *p_term = results_array;
+ FComplexe *p_theta_derivated_term = theta_derivated_results_array;
+ double *p_spherical_harmonic_Inner_coefficients_array = sphereHarmoInnerCoef;
+ double *ptr_associated_Legendre_function_Array = legendre;
+ double *start_ptr_associated_Legendre_function_Array = ptr_associated_Legendre_function_Array;
+
+
+ // Initialization of precomputed_cos_and_sin_array:
+ for(int idxm = 0 , idxmMod4 = 0; idxm <= FMB_Info_P ; ++idxm, ++idxmMod4){
+ if(idxmMod4 == 4) idxmMod4 = 0;
+ const double angle = idxm*inSphere.phi + PiArrayInner[idxmMod4];
+ (cosSin + idxm)->setReal(FMath::Sin(angle + M_PI_2));
+ (cosSin + idxm)->setImag(FMath::Sin(angle));
+ }
+
+ // Initialization of associated_Legendre_function_Array:
+ legendreFunction( inSphere.cosTheta, inSphere.sinTheta, this->legendre);
+
+ // r^l
+ double r_l = 1.0;
+ for (int l = 0 ; l <= FMB_Info_P ; ++l, r_l *= inSphere.r){
+ double magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (*ptr_associated_Legendre_function_Array);
+ // m<l:
+ int m = 0;
+ for(; m < l ; ++m, ++p_term, ++p_theta_derivated_term, ++p_spherical_harmonic_Inner_coefficients_array, ++ptr_associated_Legendre_function_Array){
+ // Computation of Inner_l^m(r, theta, phi):
+ p_term->setReal( magnitude * (cosSin + m)->getReal());
+ p_term->setImag( magnitude * (cosSin + m)->getImag());
+ // Computation of {\partial Inner_l^m(r, theta, phi)}/{\partial theta}:
+ magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * ((l*inSphere.cosTheta*(*ptr_associated_Legendre_function_Array)
+ - (l+m)*(*(start_ptr_associated_Legendre_function_Array + expansion_Redirection_array_for_j[l-1] + m))) / inSphere.sinTheta);
+ p_theta_derivated_term->setReal(magnitude * (cosSin + m)->getReal());
+ p_theta_derivated_term->setImag(magnitude * (cosSin + m)->getImag());
+ }
+
+ // m=l:
+ // Computation of Inner_m^m(r, theta, phi):
+ magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (*ptr_associated_Legendre_function_Array);
+ p_term->setReal(magnitude * (cosSin + m)->getReal());
+ p_term->setImag(magnitude * (cosSin + m)->getImag());
+ // Computation of {\partial Inner_m^m(r, theta, phi)}/{\partial theta}:
+ magnitude = (*p_spherical_harmonic_Inner_coefficients_array) * r_l * (m * inSphere.cosTheta * (*ptr_associated_Legendre_function_Array) / inSphere.sinTheta);
+ p_theta_derivated_term->setReal(magnitude * (cosSin + m)->getReal());
+ p_theta_derivated_term->setImag(magnitude * (cosSin + m)->getImag());
+
+ ++p_term;
+ ++p_theta_derivated_term;
+ ++p_spherical_harmonic_Inner_coefficients_array;
+ ++ptr_associated_Legendre_function_Array;
+ }
+ }
+
+ //[fmb] ff_matrix_Convert_exp_2_transfer_M2L_matrix
+ /*[Blas] void expffmatrix(const FComplexe* transfer_exp, FComplexe* const ff_transfer_matrix){
+ FComplexe *p_ff_transfer_matrix = ff_transfer_matrix;
+
+ for( int M = 0; M <= this->FMB_Info_P ; ++M ){
+ for(int m = 0; m <= M ; ++m ){
+ for(int N = 0; N <= FMB_Info_P ; ++N ){
+ for(int n = 0; n <= 2 * N ; ++n , ++p_ff_transfer_matrix ){
+ const int k = N - n - m;
+ if(k < 0){
+ const int pow_of_minus_1_k = ((k%2) ? -1 : 1);
+ p_ff_transfer_matrix->setReal(pow_of_minus_1_k * transfer_exp[ expansion_Redirection_array_for_j[M + N] - k].getReal());
+ p_ff_transfer_matrix->setImag((-pow_of_minus_1_k) * transfer_exp[ expansion_Redirection_array_for_j[M + N] - k].getImag());
+ }
+ else{
+ int temp = expansion_Redirection_array_for_j[M + N];
+ *p_ff_transfer_matrix = transfer_exp[ expansion_Redirection_array_for_j[M + N] + k];
+ }
+ }
+ }
+ }
+ }
+ }*/
+
+ //////////////////////////////////////////////////////////////////
+ // Precompute
+ //////////////////////////////////////////////////////////////////
+
+ // transfer_M2M_Precompute_all_levels
+ // transfer_L2L_Precompute_all_levels
+ void precomputeM2M(){
+ double treeWidthAtLevel = this->treeWidthAtRoot/2;
+
+ for(int idxLevel = 0 ; idxLevel < this->treeHeight ; ++idxLevel ){
+ const F3DPosition father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel);
+ treeWidthAtLevel /= 2;
+
+ //std::cout << "[precomputeM2M]treeWidthAtLevel=" << treeWidthAtLevel << "\n";
+ //printf("\tidxLevel=%d\tFather.x=%f\tFather.y=%f\tFather.z=%f\n",idxLevel,father.getX(),father.getY(),father.getZ());
+
+ for(int idxChild = 0 ; idxChild < 8 ; ++idxChild ){
+ FTreeCoordinate childBox;
+ childBox.setPositionFromMorton(idxChild,1);
+
+ const F3DPosition M2MVector (
+ father.getX() - (treeWidthAtLevel * (1 + (childBox.getX() * 2))),
+ father.getY() - (treeWidthAtLevel * (1 + (childBox.getY() * 2))),
+ father.getZ() - (treeWidthAtLevel * (1 + (childBox.getZ() * 2)))
+ );
+ Spherical sphericalM2M;
+ positionToSphere(M2MVector,&sphericalM2M);
+ harmonicInner(sphericalM2M,this->transitionM2M[idxLevel][idxChild]);
+
+ const F3DPosition L2LVector (
+ (treeWidthAtLevel * (1 + (childBox.getX() * 2))) - father.getX(),
+ (treeWidthAtLevel * (1 + (childBox.getY() * 2))) - father.getY(),
+ (treeWidthAtLevel * (1 + (childBox.getZ() * 2))) - father.getZ()
+ );
+ Spherical sphericalL2L;
+ positionToSphere(L2LVector,&sphericalL2L);
+ harmonicInner(sphericalL2L,this->transitionL2L[idxLevel][idxChild]);
+
+ //std::cout << "[M2M_vector]" << idxLevel << "/" << idxChild << " = " << M2MVector.getX() << "/" << M2MVector.getY() << "/" << M2MVector.getZ() << "\n";
+ //std::cout << "[M2M_vectorSpherical]" << idxLevel << "/" << idxChild << " = " << sphericalM2M.r << "/" << sphericalM2M.cosTheta << "/" << sphericalM2M.sinTheta << "/" << sphericalM2M.phi << "\n";
+ //for(int idxExpSize = 0 ; idxExpSize < FMB_Info_exp_size ; ++idxExpSize){
+ //std::cout << "transitionL2L[" << idxLevel << "][" << idxChild << "][" << idxExpSize << "]=" << this->transitionL2L[idxLevel][idxChild][idxExpSize].getReal()<<"/"<<this->transitionL2L[idxLevel][idxChild][idxExpSize].getImag()<< "\n";
+ //std::cout << "transitionM2M[" << idxLevel << "][" << idxChild << "][" << idxExpSize << "]=" << this->transitionM2M[idxLevel][idxChild][idxExpSize].getReal()<<"/"<<this->transitionM2M[idxLevel][idxChild][idxExpSize].getImag()<< "\n";
+ //}
+ }
+ }
+
+ }
+
+
+ // transfer_M2L_Precompute_all_levels
+ void precomputeM2L(){
+ //[Blas] FComplexe tempComplexe[FMB_Info_M2L_exp_size];
+ //printf("FMB_Info.M2L_exp_size = %d\n",FMB_Info_M2L_exp_size);
+
+ double treeWidthAtLevel = this->treeWidthAtRoot;
+ for(int idxLevel = 0 ; idxLevel < this->treeHeight ; ++idxLevel ){
+ //printf("level = %d \t width = %lf\n",idxLevel,treeWidthAtLevel);
+ for( int idxd1 = 0; idxd1 < this->size1Dim ; ++idxd1 ){
+
+ for( int idxd2 = 0; idxd2 < this->size1Dim ; ++idxd2 ){
+
+ for( int idxd3 = 0; idxd3 < this->size1Dim ; ++idxd3 ){
+ const long x = idxd1 - this->halphSize1Dim;
+ const long y = idxd2 - this->halphSize1Dim;
+ const long z = idxd3 - this->halphSize1Dim;
+
+ //printf("x=%ld \t y=%ld \t z=%ld\n",x,y,z);
+
+ if( ( x*x + y*y + z*z ) >= ( 3*FMB_Info_ws*FMB_Info_ws + 0.1 ) ){
+ const F3DPosition relativePos( x*treeWidthAtLevel , y*treeWidthAtLevel , z*treeWidthAtLevel );
+
+ Spherical spherical;
+ positionToSphere(relativePos,&spherical);
+ // Not blas so
+ //printf("transferM2L[%d][%d][%d][%d]\n", idxLevel, idxd1, idxd2, idxd3);
+ harmonicOuter(spherical,this->transferM2L[idxLevel][idxd1][idxd2][idxd3]);
+ //for(int idxTemp = 0 ; idxTemp < this->FMB_Info_M2L_exp_size ; ++idxTemp){
+ // printf("transferM2L[%d][%d][%d][%d][%d]=%f/%f\n", idxLevel, idxd1, idxd2, idxd3, idxTemp, this->transferM2L[idxLevel][idxd1][idxd2][idxd3][idxTemp].getReal(),this->transferM2L[idxLevel][idxd1][idxd2][idxd3][idxTemp].getImag());
+ //}
+
+ //[Blas] harmonicOuter(spherical,tempComplexe);
+ //[Blas] ff_matrix_Convert_exp_2_transfer_M2L_matrix
+ //[Blas] expffmatrix( tempComplexe , this->transferM2L[idxLevel][idxd1][idxd2][idxd3] );
+ }
+
+ }
+ }
+ }
+ treeWidthAtLevel /= 2;
+ }
+ }
+
+ void buildPrecompute(){
+ expansion_Redirection_array_for_j_Initialize();
+ sphericalHarmonicInitialize();
+ transferAllocate();
+
+ precomputeM2M();
+ precomputeM2L();
+ }
+
+public:
+ FFmbKernels(const int inTreeHeight, const double inTreeWidth)
+ : treeHeight(inTreeHeight), treeWidthAtRoot(inTreeWidth),
+ transitionM2M(0), transitionL2L(0), potential_sum(0){
+ buildPrecompute();
+ }
+
+ /** Default destructor */
+ virtual ~FFmbKernels(){
+ transferDeallocate();
+ }
+
+ F3DPosition getForcesSum() const {
+ return force_sum;
+ }
+
+ double getPotential() const{
+ return potential_sum;
+ }
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Upward
+ /////////////////////////////////////////////////////////////////////////////////
+
+ /** OK!
+ * expansion_P2M_add
+ * Multipole expansion with m charges q_i in Q_i=(rho_i, alpha_i, beta_i)
+ *whose relative coordinates according to *p_center are:
+ *Q_i - *p_center = (rho'_i, alpha'_i, beta'_i);
+ *
+ *For j=0..P, k=-j..j, we have:
+ *
+ *M_j^k = (-1)^j { sum{i=1..m} q_i Inner_j^k(rho'_i, alpha'_i, beta'_i) }
+ *
+ *However the extern loop is over the bodies (i=1..m) in our code and as an
+ *intern loop we have: j=0..P, k=-j..j
+ *
+ *and the potential is then given by:
+ *
+ * Phi(x) = sum_{n=0}^{+} sum_{m=-n}^{n} M_n^m O_n^{-m} (x - *p_center)
+ *
+ */
+ void P2M(CellClass* const inPole, FList<ParticuleClass*>* const inParticules) {
+
+ for(typename FList<ParticuleClass*>::BasicIterator iterParticule(*inParticules);
+ iterParticule.isValide() ; iterParticule.progress()){
+
+ Spherical spherical;
+ positionToSphere(iterParticule.value()->getPosition() - inPole->getPosition(), &spherical);
+
+ //std::cout << "Working on part " << iterParticule.value()->getValue() << "\n";
+ F3DPosition tempPos = iterParticule.value()->getPosition() - inPole->getPosition();
+ //ok printf("\tpos_rel.x=%f\tpos_rel.y=%f\tpos_rel.z=%f\n",tempPos.getX(),tempPos.getY(),tempPos.getZ());
+ //ok printf("\tp_center.x=%f\tp_center.y=%f\tp_center.z=%f\n",inPole->getPosition().getX(),inPole->getPosition().getY(),inPole->getPosition().getZ());
+ //ok printf("\tbody.x=%f\tbody.y=%f\tbody.z=%f\n",iterParticule.value()->getPosition().getX(),iterParticule.value()->getPosition().getY(),iterParticule.value()->getPosition().getZ());
+
+ harmonicInner(spherical,current_thread_Y);
+
+ //ok printf("\tr=%f\tcos_theta=%f\tsin_theta=%f\tphi=%f\n",spherical.r,spherical.cosTheta,spherical.sinTheta,spherical.phi);
+
+ FComplexe* p_exp_term = inPole->getMultipole();
+ FComplexe* p_Y_term = current_thread_Y;
+ double pow_of_minus_1_j = 1.0;//(-1)^j
+ const double valueParticule = iterParticule.value()->getValue();
+
+ for(int j = 0 ; j <= FMB_Info_P ; ++j, pow_of_minus_1_j = -pow_of_minus_1_j ){
+ for(int k = 0 ; k <= j ; ++k, ++p_Y_term, ++p_exp_term){
+ p_Y_term->mulRealAndImag( valueParticule * pow_of_minus_1_j );
+ (*p_exp_term) += (*p_Y_term);
+ //printf("\tj=%d\tk=%d\tp_exp_term.real=%f\tp_exp_term.imag=%f\tp_Y_term.real=%f\tp_Y_term.imag=%f\tpow_of_minus_1_j=%f\n", j,k,(*p_exp_term).getReal(),(*p_exp_term).getImag(),(*p_Y_term).getReal(),(*p_Y_term).getImag(),pow_of_minus_1_j);
+ }
+ }
+ }
+ }
+
+ /**
+ *-----------------------------------
+ *octree_Upward_pass_internal_cell
+ *expansion_M2M_add
+ *-----------------------------------
+ *We compute the translation of multipole_exp_src from *p_center_of_exp_src to
+ *p_center_of_exp_target, and add the result to multipole_exp_target.
+ *
+ * O_n^l (with n=0..P, l=-n..n) being the former multipole expansion terms
+ * (whose center is *p_center_of_multipole_exp_src) we have for the new multipole
+ * expansion terms (whose center is *p_center_of_multipole_exp_target):
+
+ * M_j^k = sum{n=0..j}
+ * sum{l=-n..n, |k-l|<=j-n}
+ * O_n^l Inner_{j-n}^{k-l}(rho, alpha, beta)
+ *
+ * where (rho, alpha, beta) are the spherical coordinates of the vector :
+ * p_center_of_multipole_exp_target - *p_center_of_multipole_exp_src
+ *
+ * Warning: if j-n < |k-l| we do nothing.
+ */
+ void M2M(CellClass* const inPole, CellClass** const inChild, const int inLevel) {
+
+ // We do NOT have: for(l=n-j+k; l<=j-n+k ;++l){} <=> for(l=-n; l<=n ;++l){if (j-n >= abs(k-l)){}}
+ // But we have: for(k=MAX(0,n-j+l); k<=j-n+l; ++k){} <=> for(k=0; k<=j; ++k){if (j-n >= abs(k-l)){}}
+ // (This is not the same as in L2L since the range of values of k is not the same, compared to "n-j" or "j-n".)
+ // Therefore the loops over n and l are the outmost ones and
+ // we invert the loop over j with the summation with n:
+ // for{j=0..P} sum_{n=0}^j <-> sum_{n=0}^P for{j=n..P}
+ FComplexe* const multipole_exp_target = inPole->getMultipole();
+
+ for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
+ if(!inChild[idxChild]) continue;
+ //printf("\tChild %d\n",idxChild);
+
+ const FComplexe* const multipole_exp_src = inChild[idxChild]->getMultipole();
+
+ const FComplexe* const M2M_transfer = transitionM2M[inLevel][idxChild];
+
+ for(int n = 0 ; n <= FMB_Info_P ; ++n ){
+ // l<0 // (-1)^l
+ double pow_of_minus_1_for_l = ( n % 2 ? -1 : 1);
+
+ // O_n^l : here points on the source multipole expansion term of degree n and order |l|
+ const FComplexe* p_src_exp_term = multipole_exp_src + expansion_Redirection_array_for_j[n]+n;
+ //printf("\t[p_src_exp_term] expansion_Redirection_array_for_j[n]=%d\tn=%d\n",expansion_Redirection_array_for_j[n],n);
+
+ int l = -n;
+ for(; l<0 ; ++l, --p_src_exp_term, pow_of_minus_1_for_l = -pow_of_minus_1_for_l){
+
+ for(int j = n ; j<= FMB_Info_P ; ++j ){
+ // M_j^k
+ FComplexe *p_target_exp_term = multipole_exp_target + expansion_Redirection_array_for_j[j];
+ //printf("\t[p_target_exp_term] expansion_Redirection_array_for_j[j]=%d\n",expansion_Redirection_array_for_j[j]);
+ // Inner_{j-n}^{k-l} : here points on the M2M transfer function/expansion term of degree n-j and order |k-l|
+ const FComplexe *p_Inner_term= M2M_transfer + expansion_Redirection_array_for_j[j-n]-l /* k==0 */;
+ //printf("\t[p_Inner_term] expansion_Redirection_array_for_j[j-n]=%d\tl=%d\n",expansion_Redirection_array_for_j[j-n],-l);
+
+ // since n-j+l<0
+ for(int k=0 ; k <= (j-n+l) ; ++k, ++p_target_exp_term, ++p_Inner_term){ // l<0 && k>=0 => k-l>0
+ p_target_exp_term->incReal( pow_of_minus_1_for_l *
+ ((p_src_exp_term->getReal() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getImag() * p_Inner_term->getImag())));
+ p_target_exp_term->incImag( pow_of_minus_1_for_l *
+ ((p_src_exp_term->getReal() * p_Inner_term->getImag()) -
+ (p_src_exp_term->getImag() * p_Inner_term->getReal())));
+
+ //printf("\tp_src_exp_term->getReal()=%f\tp_src_exp_term->getImag()=%f\n", p_src_exp_term->getReal(),p_src_exp_term->getImag());
+ //printf("\tp_Inner_term->getReal()=%f\tp_Inner_term->getImag()=%f\n", p_Inner_term->getReal(),p_Inner_term->getImag());
+ //printf("\tn[%d]l[%d]j[%d]k[%d] = %f / %f\n",n,l,j,k,p_target_exp_term->getReal(),p_target_exp_term->getImag());
+
+ } // for k
+ } // for j
+ } // for l
+
+ // l>=0
+ for(; l <= n ; ++l, ++p_src_exp_term, pow_of_minus_1_for_l = -pow_of_minus_1_for_l){
+
+ for( int j=n ; j <= FMB_Info_P ; ++j ){
+ // (-1)^k
+ double pow_of_minus_1_for_k = ( FMath::Max(0,n-j+l) %2 ? -1 : 1 );
+ // M_j^k
+ FComplexe *p_target_exp_term = multipole_exp_target + expansion_Redirection_array_for_j[j] + FMath::Max(0,n-j+l);
+ // Inner_{j-n}^{k-l} : here points on the M2M transfer function/expansion term of degree n-j and order |k-l|
+ const FComplexe *p_Inner_term = M2M_transfer + expansion_Redirection_array_for_j[j-n] + l - FMath::Max(0,n-j+l);// -(k-l)
+
+ int k = FMath::Max(0,n-j+l);
+ for(; k <= (j-n+l) && (k-l) < 0 ; ++k, ++p_target_exp_term, --p_Inner_term, pow_of_minus_1_for_k = -pow_of_minus_1_for_k){ /* l>=0 && k-l<0 */
+ p_target_exp_term->incReal( pow_of_minus_1_for_k * pow_of_minus_1_for_l *
+ ((p_src_exp_term->getReal() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getImag() * p_Inner_term->getImag())));
+ p_target_exp_term->incImag(pow_of_minus_1_for_k * pow_of_minus_1_for_l *
+ ((p_src_exp_term->getImag() * p_Inner_term->getReal()) -
+ (p_src_exp_term->getReal() * p_Inner_term->getImag())));
+ //printf("\tn[%d]l[%d]j[%d]k[%d] = %f / %f\n",n,l,j,k,p_target_exp_term->getReal(),p_target_exp_term->getImag());
+
+ } // for k
+
+ for(; k <= (j - n + l) ; ++k, ++p_target_exp_term, ++p_Inner_term){ // l>=0 && k-l>=0
+ p_target_exp_term->incReal(
+ (p_src_exp_term->getReal() * p_Inner_term->getReal()) -
+ (p_src_exp_term->getImag() * p_Inner_term->getImag()));
+ p_target_exp_term->incImag(
+ (p_src_exp_term->getImag() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getReal() * p_Inner_term->getImag()));
+ //printf("\tn[%d]l[%d]j[%d]k[%d] = %f / %f\n",n,l,j,k,p_target_exp_term->getReal(),p_target_exp_term->getImag());
+
+ } // for k
+ } // for j
+ } // for l
+ } // for n
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // M2L
+ /////////////////////////////////////////////////////////////////////////////////
+ /**
+ *------------------
+ * expansion_M2L_add
+ *-------------------
+ *We compute the conversion of multipole_exp_src in *p_center_of_exp_src to
+ *a local expansion in *p_center_of_exp_target, and add the result to local_exp_target.
+ *
+ *O_n^l (with n=0..P, l=-n..n) being the former multipole expansion terms
+ *(whose center is *p_center_of_multipole_exp_src) we have for the new local
+ *expansion terms (whose center is *p_center_of_local_exp_target):
+ *
+ *L_j^k = sum{n=0..+}
+ *sum{l=-n..n}
+ *O_n^l Outer_{j+n}^{-k-l}(rho, alpha, beta)
+ *
+ *where (rho, alpha, beta) are the spherical coordinates of the vector :
+ *p_center_of_local_exp_src - *p_center_of_multipole_exp_target
+ *
+ *Remark: here we have always j+n >= |-k-l|
+ *
+ */
+ void M2L(CellClass* const pole, CellClass** const distantNeighbors, const int size, const int inLevel) {
+ FTreeCoordinate coordCenter;
+ coordCenter.setPositionFromMorton(pole->getMortonIndex(),inLevel);
+
+ bool FMB_Info_up_to_P_in_M2L = true;
+
+ for(int idxSize = 0 ; idxSize < size ; ++idxSize){
+ FTreeCoordinate coordNeighbors;
+ coordNeighbors.setPositionFromMorton(distantNeighbors[idxSize]->getMortonIndex(),inLevel);
+
+ //printf("\tidxSize = %d\tleve = %d\tMorton = %lld\n",idxSize,inLevel,distantNeighbors[idxSize]->getMortonIndex());
+
+ const FComplexe* const M2L_transfer = transferM2L[inLevel]
+ [6-(coordNeighbors.getX()+halphSize1Dim-coordCenter.getX())]
+ [6-(coordNeighbors.getY()+halphSize1Dim-coordCenter.getY())]
+ [6-(coordNeighbors.getZ()+halphSize1Dim-coordCenter.getZ())];
+ /*printf("level = %d\tx=%ld\ty=%ld\tz=%ld\n", inLevel,
+ 6-(coordNeighbors.getX()+halphSize1Dim-coordCenter.getX()),
+ 6-(coordNeighbors.getY()+halphSize1Dim-coordCenter.getY()),
+ 6-(coordNeighbors.getZ()+halphSize1Dim-coordCenter.getZ()));*/
+ /*printf("M2L_transfer[0]= %f/%f\n",M2L_transfer->getReal(),M2L_transfer->getImag());
+ printf("M2L_transfer[1]= %f/%f\n",M2L_transfer[1].getReal(),M2L_transfer[1].getImag());
+ printf("M2L_transfer[2]= %f/%f\n",M2L_transfer[2].getReal(),M2L_transfer[2].getImag());*/
+ const FComplexe* const multipole_exp_src = distantNeighbors[idxSize]->getMultipole();
+ // L_j^k
+ FComplexe* p_target_exp_term = pole->getLocal();
+ int start_for_j = 0;
+ //#if defined (_FORCES_) && !defined(_ENERGY_)
+ // See FMB.c:
+ if (FMB_Info_up_to_P_in_M2L){
+ start_for_j = 1;
+ ++p_target_exp_term;
+ }
+ //#endif // #if defined (_FORCES_) && !defined(_ENERGY_)
+
+ // HPMSTART(51, "M2L computation (loops)");
+ for (int j = start_for_j ; j <= FMB_Info_P ; ++j){
+
+ int stop_for_n = FMB_Info_P;
+ if (FMB_Info_up_to_P_in_M2L){
+ stop_for_n = FMB_Info_P - j;
+ }
+
+ // (-1)^k
+ double pow_of_minus_1_for_k = 1.0;
+ for (int k = 0 ; k <= j ; ++k, pow_of_minus_1_for_k = -pow_of_minus_1_for_k, ++p_target_exp_term){
+
+ // (-1)^n
+ double pow_of_minus_1_for_n = 1.0;
+ for (int n = 0 ; n <= stop_for_n ; ++n, pow_of_minus_1_for_n = -pow_of_minus_1_for_n){
+
+ // O_n^l : here points on the source multipole expansion term of degree n and order |l|
+ const FComplexe *p_src_exp_term = multipole_exp_src + expansion_Redirection_array_for_j[n] + n;
+ // Outer_{j+n}^{-k-l} : here points on the M2L transfer function/expansion term of degree j+n and order |-k-l|
+ const FComplexe *p_Outer_term = M2L_transfer + expansion_Redirection_array_for_j[n+j] + k+n;
+ //printf("expansion_Get_p_term(M2L_transfer, n+j, k+n)-M2L_transfer = %d \t(n=%d)\n",
+ // p_Outer_term-M2L_transfer , n );
+ //printf("TRUC = %d\n",expansion_Redirection_array_for_j[n+j] + k+n);
+ double pow_of_minus_1_for_l = pow_of_minus_1_for_n; // (-1)^l
+ // We start with l=n (and not l=-n) so that we always set p_Outer_term to a correct value in the first loop.
+ int l=n;
+ for ( ; l>0 ; --l, pow_of_minus_1_for_l = -pow_of_minus_1_for_l, --p_src_exp_term, --p_Outer_term){ // we have -k-l<0 and l>0
+ p_target_exp_term->incReal( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
+ ((p_src_exp_term->getReal() * p_Outer_term->getReal()) +
+ (p_src_exp_term->getImag() * p_Outer_term->getImag())));
+ p_target_exp_term->incImag( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
+ ((p_src_exp_term->getImag() * p_Outer_term->getReal()) -
+ (p_src_exp_term->getReal() * p_Outer_term->getImag())));
+ //printf("\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ // p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ //printf("\t p_src_exp_term->real = %lf \t p_src_exp_term->imag = %lf \n",
+ // p_src_exp_term->getReal(),p_src_exp_term->getImag());
+ //printf("p_Outer_term-M2L_transfer = %d\n",
+ // p_Outer_term-M2L_transfer);
+ //printf("\t p_Outer_term->real = %f \t p_Outer_term->imag = %f \n",
+ // p_Outer_term->getReal(),p_Outer_term->getImag());
+ }
+
+ for (; l>=-n && -k-l<0 ; --l, pow_of_minus_1_for_l = -pow_of_minus_1_for_l, ++p_src_exp_term, --p_Outer_term){ // we have -k-l<0 and l<=0
+ p_target_exp_term->incReal( pow_of_minus_1_for_k *
+ ((p_src_exp_term->getReal() * p_Outer_term->getReal()) -
+ (p_src_exp_term->getImag() * p_Outer_term->getImag())));
+ p_target_exp_term->decImag( pow_of_minus_1_for_k *
+ ((p_src_exp_term->getImag() * p_Outer_term->getReal()) +
+ (p_src_exp_term->getReal() * p_Outer_term->getImag())));
+ }
+
+ for (; l>=-n; --l, pow_of_minus_1_for_l = -pow_of_minus_1_for_l, ++p_src_exp_term, ++p_Outer_term){ // we have -k-l>=0 and l<=0
+ p_target_exp_term->incReal( pow_of_minus_1_for_l *
+ ((p_src_exp_term->getReal() * p_Outer_term->getReal()) +
+ (p_src_exp_term->getImag() * p_Outer_term->getImag())));
+ p_target_exp_term->incImag( pow_of_minus_1_for_l *
+ ((p_src_exp_term->getReal() * p_Outer_term->getImag()) -
+ (p_src_exp_term->getImag() * p_Outer_term->getReal())));
+ }
+ //printf("\tj=%d\tk=%d\tn=%d\tl=%d\n",j,k,n,l);
+ //printf("\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ // p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ }
+ }
+ }
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Downard
+ /////////////////////////////////////////////////////////////////////////////////
+
+ /**
+ *We compute the shift of local_exp_src from *p_center_of_exp_src to
+ *p_center_of_exp_target, and set the result to local_exp_target.
+ *
+ *O_n^l (with n=0..P, l=-n..n) being the former local expansion terms
+ *(whose center is *p_center_of_exp_src) we have for the new local
+ *expansion terms (whose center is *p_center_of_exp_target):
+ *
+ *L_j^k = sum{n=j..P}
+ *sum{l=-n..n}
+ *O_n^l Inner_{n-j}^{l-k}(rho, alpha, beta)
+ *
+ *where (rho, alpha, beta) are the spherical coordinates of the vector :
+ *p_center_of_exp_target - *p_center_of_exp_src
+ *
+ *Warning: if |l-k| > n-j, we do nothing.
+ */
+ void L2L(CellClass* const pole, CellClass** const child, const int inLevel) {
+
+ for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
+ // if no child at this position
+ if(!child[idxChild]) continue;
+
+ const FComplexe* const L2L_tranfer = transitionL2L[inLevel][idxChild];
+ const FComplexe* const local_exp_src = pole->getLocal();
+ FComplexe* const local_exp_target = child[idxChild]->getLocal();
+
+ /*printf("Level %d\n", inLevel);
+ printf("Father morton %lld\n", pole->getMortonIndex());
+ printf("Child morton %lld\n", child[idxChild]->getMortonIndex());*/
+
+ // L_j^k
+ FComplexe* p_target_exp_term = local_exp_target;
+ for (int j=0 ; j<= FMB_Info_P ; ++j){
+ // (-1)^k
+ double pow_of_minus_1_for_k = 1.0;
+ for (int k=0 ; k <= j ; ++k, pow_of_minus_1_for_k = -pow_of_minus_1_for_k, ++p_target_exp_term){
+ for (int n=j; n<=FMB_Info_P;++n){
+ // O_n^l : here points on the source multipole expansion term of degree n and order |l|
+ const FComplexe* p_src_exp_term = local_exp_src + expansion_Redirection_array_for_j[n] + n-j+k;
+ //printf("expansion_Redirection_array_for_j[n] + n-j+k %d\n", expansion_Redirection_array_for_j[n] + n-j+k);
+ int l = n-j+k;
+ // Inner_{n-j}^{l-k} : here points on the L2L transfer function/expansion term of degree n-j and order |l-k|
+ const FComplexe* p_Inner_term = L2L_tranfer + expansion_Redirection_array_for_j[n-j] + l-k;
+
+ //printf("1\n");
+ for ( ; l-k>0; --l, --p_src_exp_term, --p_Inner_term){ /* l>0 && l-k>0 */
+ p_target_exp_term->incReal( (p_src_exp_term->getReal() * p_Inner_term->getReal()) -
+ (p_src_exp_term->getImag() * p_Inner_term->getImag()));
+ p_target_exp_term->incImag( (p_src_exp_term->getImag() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getReal() * p_Inner_term->getImag()));
+ /*printf("\t p_src_exp_term->real = %lf \t p_src_exp_term->imag = %lf \n",
+ p_src_exp_term->getReal(),p_src_exp_term->getImag());
+ printf("\t p_Inner_term->real = %lf \t p_Inner_term->imag = %lf \n",
+ p_Inner_term->getReal(),p_Inner_term->getImag());
+ printf("\t\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ printf("\tp_target_exp_term = %d\n",p_target_exp_term-local_exp_target);*/
+ }
+
+ //printf("2\n");
+ // (-1)^l
+ double pow_of_minus_1_for_l = ((l%2) ? -1 : 1);
+ for (; l>0 && l>=j-n+k; --l, pow_of_minus_1_for_l = -pow_of_minus_1_for_l, --p_src_exp_term, ++p_Inner_term){ /* l>0 && l-k<=0 */
+ p_target_exp_term->incReal( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
+ ((p_src_exp_term->getReal() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getImag() * p_Inner_term->getImag())));
+ p_target_exp_term->incImag( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
+ ((p_src_exp_term->getImag() * p_Inner_term->getReal()) -
+ (p_src_exp_term->getReal() * p_Inner_term->getImag())));
+ /*printf("\t p_src_exp_term->real = %lf \t p_src_exp_term->imag = %lf \n",
+ p_src_exp_term->getReal(),p_src_exp_term->getImag());
+ printf("\t p_Inner_term->real = %lf \t p_Inner_term->imag = %lf \n",
+ p_Inner_term->getReal(),p_Inner_term->getImag());
+ printf("\t\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ printf("\tp_target_exp_term = %d\n",p_target_exp_term-local_exp_target);*/
+ }
+
+ //printf("3\n");
+ // l<=0 && l-k<=0
+ for (; l>=j-n+k; --l, ++p_src_exp_term, ++p_Inner_term){
+ p_target_exp_term->incReal( pow_of_minus_1_for_k *
+ ((p_src_exp_term->getReal() * p_Inner_term->getReal()) -
+ (p_src_exp_term->getImag() * p_Inner_term->getImag())));
+ p_target_exp_term->decImag( pow_of_minus_1_for_k *
+ ((p_src_exp_term->getImag() * p_Inner_term->getReal()) +
+ (p_src_exp_term->getReal() * p_Inner_term->getImag())));
+ /*printf("\t p_src_exp_term->real = %lf \t p_src_exp_term->imag = %lf \n",
+ p_src_exp_term->getReal(),p_src_exp_term->getImag());
+ printf("\t p_Inner_term->real = %lf \t p_Inner_term->imag = %lf \n",
+ p_Inner_term->getReal(),p_Inner_term->getImag());
+ printf("\t\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ printf("\tj=%d\tk=%d\tn=%d\tl=%d\tpow_of_minus_1_for_k=%f\n",j,k,n,l,pow_of_minus_1_for_k);
+ printf("\tp_target_exp_term = %d\n",p_target_exp_term-local_exp_target);*/
+ }
+ //printf("\tj=%d\tk=%d\tn=%d\tl=%d\n",j,k,n,l);
+ //printf("\t\t p_target_exp_term->real = %lf \t p_target_exp_term->imag = %lf \n",
+ // p_target_exp_term->getReal(),p_target_exp_term->getImag());
+ }
+ }
+ }
+ }
+ }
+
+ /** bodies_L2P
+ * expansion_L2P_add_to_force_vector
+ */
+ void L2P(CellClass* const local, FList<ParticuleClass*>* const particules){
+ typename FList<ParticuleClass*>::BasicIterator iterTarget(*particules);
+ while( iterTarget.isValide() ){
+
+ F3DPosition force_vector_in_local_base;
+ Spherical spherical;
+
+ positionToSphere( iterTarget.value()->getPosition() - local->getPosition(), &spherical );
+
+ harmonicInnerThetaDerivated( spherical, current_thread_Y, current_thread_Y_theta_derivated);
+
+ // The maximum degree used here will be P.
+ FComplexe* p_Y_term = current_thread_Y+1;
+ FComplexe* p_Y_theta_derivated_term = current_thread_Y_theta_derivated+1;
+ FComplexe* p_local_exp_term = local->getLocal()+1;
+ for (int j = 1 ; j <= FMB_Info_P ; ++j ){
+ FComplexe exp_term_aux;
+
+ // k=0:
+ // F_r:
+ exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
+ exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
+
+ force_vector_in_local_base.setX( force_vector_in_local_base.getX() + j * exp_term_aux.getReal());
+ // F_phi: k=0 => nothing to do for F_phi
+ // F_theta:
+ exp_term_aux.setReal( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
+ exp_term_aux.setImag( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getImag()) + (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getReal()) );
+
+ force_vector_in_local_base.setY( force_vector_in_local_base.getY() + exp_term_aux.getReal());
+
+ ++p_local_exp_term;
+ ++p_Y_term;
+ ++p_Y_theta_derivated_term;
+
+
+ // k>0:
+ for (int k=1; k<=j ;++k, ++p_local_exp_term, ++p_Y_term, ++p_Y_theta_derivated_term){
+ // F_r:
+
+ exp_term_aux.setReal( (p_Y_term->getReal() * p_local_exp_term->getReal()) - (p_Y_term->getImag() * p_local_exp_term->getImag()) );
+ exp_term_aux.setImag( (p_Y_term->getReal() * p_local_exp_term->getImag()) + (p_Y_term->getImag() * p_local_exp_term->getReal()) );
+
+ force_vector_in_local_base.setX(force_vector_in_local_base.getX() + 2 * j * exp_term_aux.getReal());
+ // F_phi:
+ force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() - 2 * k * exp_term_aux.getImag());
+ // F_theta:
+
+ exp_term_aux.setReal( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getReal()) - (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getImag()) );
+ exp_term_aux.setImag( (p_Y_theta_derivated_term->getReal() * p_local_exp_term->getImag()) + (p_Y_theta_derivated_term->getImag() * p_local_exp_term->getReal()) );
+
+ force_vector_in_local_base.setY(force_vector_in_local_base.getY() + 2 * exp_term_aux.getReal());
+ }
+ }
+
+ // We want: - gradient(POTENTIAL_SIGN potential).
+ // The -(- 1.0) computing is not the most efficient programming ...
+ //#define FMB_TMP_SIGN -(POTENTIAL_SIGN 1.0)
+ force_vector_in_local_base.setX( force_vector_in_local_base.getX() * (-1.0) / spherical.r);
+ force_vector_in_local_base.setY( force_vector_in_local_base.getY() * (-1.0) / spherical.r);
+ force_vector_in_local_base.setZ( force_vector_in_local_base.getZ() * (-1.0) / (spherical.r * spherical.sinTheta));
+ //#undef FMB_TMP_SIGN
+
+ const double th = FMath::ACos(spherical.cosTheta);
+ const double cos_theta = FMath::Cos(th);
+ const double cos_phi = FMath::Cos(spherical.phi);
+ const double sin_theta = FMath::Sin(th);
+ const double sin_phi = FMath::Sin(spherical.phi);
+
+ F3DPosition force_vector_tmp;
+
+ force_vector_tmp.setX(
+ cos_phi * sin_theta * force_vector_in_local_base.getX() +
+ cos_phi * cos_theta * force_vector_in_local_base.getY() +
+ (-sin_phi) * force_vector_in_local_base.getZ());
+
+ force_vector_tmp.setY(
+ sin_phi * sin_theta * force_vector_in_local_base.getX() +
+ sin_phi * cos_theta * force_vector_in_local_base.getY() +
+ cos_phi * force_vector_in_local_base.getZ());
+
+ force_vector_tmp.setZ(
+ cos_theta * force_vector_in_local_base.getX() +
+ (-sin_theta) * force_vector_in_local_base.getY());
+
+
+ //#ifndef _DIRECT_MATRIX_
+ // When _DIRECT_MATRIX_ is defined, this multiplication is done in 'leaf_Sum_near_and_far_fields()'
+ force_vector_tmp *= iterTarget.value()->getValue();
+ //#endif
+
+ iterTarget.value()->setForces( iterTarget.value()->getForces() + force_vector_tmp );
+
+ //printf("Morton %lld \t fx = %f \t fy = %f \t fz = %f \n",local->getMortonIndex(),iterTarget.value()->getForces().getX(),iterTarget.value()->getForces().getY(),iterTarget.value()->getForces().getZ());
+
+ iterTarget.progress();
+ }
+ }
+
+
+
+
+
+
+ /** void bodies_Compute_direct_interaction (
+ bodies_t *FMB_RESTRICT p_b_target,
+ bodies_t *FMB_RESTRICT p_b_src,
+ bool mutual
+ )
+ *
+ */
+ void P2P(FList<ParticuleClass*>* const currentBox, FList<ParticuleClass*>** directNeighbors, const int size) {
+ typename FList<ParticuleClass*>::BasicIterator iterTarget(*currentBox);
+ while( iterTarget.isValide() ){
+ for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < size ; ++idxDirectNeighbors){
+ typename FList<ParticuleClass*>::BasicIterator iterSource(*currentBox);
+ while( iterSource.isValide() ){
+ DIRECT_COMPUTATION_NO_MUTUAL_SOFT(&iterTarget.value(),
+ iterSource.value());
+ iterSource.progress();
+ }
+ }
+
+ typename FList<ParticuleClass*>::BasicIterator iterSameBox(*currentBox);
+ while( iterSameBox.isValide() ){
+ if(iterSameBox.value() != iterTarget.value()){
+ DIRECT_COMPUTATION_NO_MUTUAL_SOFT(&iterTarget.value(),
+ iterSameBox.value());
+ }
+ iterSameBox.progress();
+ }
+
+ //printf("x = %f \t y = %f \t z = %f \n",iterTarget.value()->getPosition().getX(),iterTarget.value()->getPosition().getY(),iterTarget.value()->getPosition().getZ());
+ //printf("\t fx = %f \t fy = %f \t fz = %f \n",iterTarget.value()->getForces().getX(),iterTarget.value()->getForces().getY(),iterTarget.value()->getForces().getZ());
+ //printf("\t potential = %f \n",iterTarget.value()->getPotential());
+
+ force_sum += iterTarget.value()->getForces();
+ potential_sum += iterTarget.value()->getPotential();
+
+ iterTarget.progress();
+ }
+ }
+ void DIRECT_COMPUTATION_NO_MUTUAL_SOFT(ParticuleClass** const target, const ParticuleClass* const source){
+ const double dx = (*target)->getPosition().getX() - source->getPosition().getX();
+ const double dy = (*target)->getPosition().getY() - source->getPosition().getY();
+ const double dz = (*target)->getPosition().getZ() - source->getPosition().getZ();
+
+ double inv_square_distance = 1.0/ (dx*dx + dy*dy + dz*dz + FMB_Info_eps_soft_square);
+ double inv_distance = FMath::Sqrt(inv_square_distance);
+ inv_distance *= (*target)->getValue() * source->getValue();
+ inv_square_distance *= inv_distance;
+
+ //#ifdef _FORCES_
+ (*target)->setForces(
+ (*target)->getForces().getX() + dx * inv_square_distance,
+ (*target)->getForces().getY() + dy * inv_square_distance,
+ (*target)->getForces().getZ() + dz * inv_square_distance
+ );
+ //#endif
+
+ //#ifdef _ENERGY_
+ (*target)->setPotential( inv_distance + (*target)->getPotential());
+ //#endif
+ }
+};
+
+
+template< class ParticuleClass, class CellClass>
+const double FFmbKernels<ParticuleClass,CellClass>::PiArrayInner[4] = {0, FMath::FPiDiv2, FMath::FPi, -FMath::FPiDiv2};
+
+
+template< class ParticuleClass, class CellClass>
+const double FFmbKernels<ParticuleClass,CellClass>::PiArrayOuter[4] = {0, -FMath::FPiDiv2, FMath::FPi, FMath::FPiDiv2};
+
+
+
+#endif //FFMBKERNELS_HPP
+
+// [--LICENSE--]