diff --git a/Doc/Src_tex/ScalFmm-PeriodicModel/periodicmodel.tex b/Doc/Src_tex/ScalFmm-PeriodicModel/periodicmodel.tex
index 3e2b60e6a8b4354f22683fab03cf850ce665c8c0..f83e9234b41b40fac760f8fe2ffcd0dc855fa808 100644
--- a/Doc/Src_tex/ScalFmm-PeriodicModel/periodicmodel.tex
+++ b/Doc/Src_tex/ScalFmm-PeriodicModel/periodicmodel.tex
@@ -138,7 +138,7 @@ Such approaches leads to a tree where all leaves may not be at the same level, b
\begin{figure}[h]
\centering
-\includegraphics[scale=0.6]{../Images/Octree}
+\includegraphics[scale=0.6]{Images/Octree}
\caption{$2D$ Octree}
\end{figure}
@@ -155,7 +155,7 @@ $r = Diagonal/2 = \sqrt{ 3 * Edge^2 }/2 = \sqrt{ 3 * (BoxWidth/2^{l})^2 }/2$.
\begin{figure}[h]
\centering
-\includegraphics[scale=0.6]{../Images/TopDownOperators}
+\includegraphics[scale=0.6]{Images/TopDownOperators}
\caption{Top-Down Operators}
\end{figure}
@@ -177,7 +177,7 @@ Any FMM implementation has the capacity to compute the $M2L$ between two cell of
\begin{figure}[h]
\centering
-\includegraphics[scale=0.6]{../Images/InteractionsOperators}
+\includegraphics[scale=0.6]{Images/InteractionsOperators}
\caption{Interactions Operators, P2P (Blue), M2L (Green)}
\end{figure}
@@ -201,7 +201,7 @@ In the next paragraph we supposed when speaking about interactions or neighbors
\begin{figure}[h]
\centering
-\includegraphics[scale=0.6]{../Images/Repetitions}
+\includegraphics[scale=0.6]{Images/Repetitions}
\caption{Repetitions of the simulation box create new neighbors}
\end{figure}
@@ -218,7 +218,7 @@ It is important to have in mind that without any special operation the box is al
\begin{figure}[h]
\centering
-\includegraphics[scale=0.35]{../Images/PeriodicL2}
+\includegraphics[scale=0.35]{Images/PeriodicL2}
\caption{Extension of the box with a periodicity until level $2$}
\end{figure}
@@ -231,7 +231,7 @@ Finally we can compute the downward pass from level $0$ to the leaves.
\begin{figure}[h]
\centering
-\includegraphics[scale=0.45]{../Images/PeriodicL0}
+\includegraphics[scale=0.45]{Images/PeriodicL0}
\caption{Extension of the box with a periodicity until level $0$}
\end{figure}
@@ -259,7 +259,7 @@ The have a perfect symmetric periodicity we need to add a border.
\begin{figure}[h]
\centering
-\includegraphics[scale=0.45]{../Images/PeriodicL-1}
+\includegraphics[scale=0.45]{Images/PeriodicL-1}
\caption{Extension of the box with a periodicity until level $-1$ without border}
\end{figure}
@@ -289,7 +289,7 @@ Each of this $M2M$ is computing the relation between one or two sub-cells and th
\begin{figure}[h]
\centering
-\includegraphics[scale=0.45]{../Images/Border}
+\includegraphics[scale=0.45]{Images/Border}
\caption{Extension of the box with a periodicity until level $-1$ with border}
\end{figure}
@@ -355,6 +355,52 @@ In fact, let the original FMM box width be $w$ and a tree height be $h$.
Then, if is required to apply periodicity until level $d$ above the root.
The kernel should be able to proceed usual FMM operator in a tree of height of size $h+d+5$ with a initial box width of $w*2^{d+3}$.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Results}
+\subsection{Ewald summation comparisons}
+
+The energy computed by molecular dynamics codes is given by
+$$
+U = \frac{1}{4 \pi\epsilon_0}\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}}
+$$
+and the force on atom $x_i$
+$$
+f(x_i) = \frac{1}{4 \pi\epsilon_0}\sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
+$$
+\subsubsection{DL\_Poly comparisons}
+DL\_POLY\_2 uses the following internal molecular units \\
+
+\begin{tabular}{|l|c|c|l|}
+\hline
+type & Expression & Numerical value & \\
+time &$t_0$ & $10^{-12}\;s$& picosecond\\
+length & $l_0$ &$10^{-10}\; m $& $\AA$ Angström\\
+mass & $m_0$ & $ 1.6605402 \; 10^{−27}\; kg $& atomic mass unit\\
+charge & $q_0$ & $1.60217733 \; 10^{−19} Coulombs$& unit of proton charge\\
+energy & $E_0 = m_0(l_0/t_0)^2$&$1.6605402 \; 10^{−23}\; Joules $ & $10\; J\, mol^{-1}$\\
+\hline
+\end{tabular}
+
+ In internal variables the energy writes
+$$
+U = \frac{q_0^2}{4 \pi\epsilon_0 l_0}\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}} = C_{energy}\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}}
+$$
+and the forces write
+$$
+f(x_i) = -\frac{q_0}{4 \pi\epsilon_0 l_0^2}\sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
+ = -C_{force}\sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
+$$
+
+The Energy conversion factor is $\gamma_0 = \frac{q_0^2}{4 \pi\epsilon_0 l_0}/E_0 = 138935.4835$. The energy unit is in Joules and if you want $kcal mol^{-1}$ unit the the factor becomes $\gamma_0/418.400$.
+\begin{tabular}{|l|c|c|c|c|c|}
+\hline
+ & \multicolumn{2}{|c|}{Energy} & \multicolumn{2}{c|}{Force} \\
+units& $kcal\, mol^{-1}$& & &\\
+\hline
+$C_{energy}$& & & &\\
+$C_{force}$& & & &\\
+\hline
+\end{tabular}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
A new approach has been proposed to compute a periodic FMM.
diff --git a/Src/Components/FAbstractParticleContainer.hpp b/Src/Components/FAbstractParticleContainer.hpp
index 14c8e5f19ef415808dd07324575747d11e1ca367..cde537f5b800252ec77c1d60b578068b7b033746 100755
--- a/Src/Components/FAbstractParticleContainer.hpp
+++ b/Src/Components/FAbstractParticleContainer.hpp
@@ -30,7 +30,7 @@ class FPoint;
* This class define the method that every particle container
* has to implement.
*
-* @warning Inherite from this class when implement a specific particle type
+* @warning Inherit from this class when implement a specific particle type
*/
class FAbstractParticleContainer {
public:
@@ -40,8 +40,8 @@ public:
/**
* This method should be inherited (or your leaf will do nothing)
- * the point is coming from the tree and is fallowed by what let the leaf
- * pass throught its push method.
+ * the point is coming from the tree and is followed by what let the leaf
+ * pass through its push method.
*/
template<typename... Args>
void push(const FPoint& /*inParticlePosition*/, Args ... /*args*/){
diff --git a/Src/Components/FBasicParticleContainer.hpp b/Src/Components/FBasicParticleContainer.hpp
index 1509eaa9cc28758331485e9b6333c87f9892273f..c4a8b2bd9da49ab4ba1ef2aa901ca9c68c808a29 100755
--- a/Src/Components/FBasicParticleContainer.hpp
+++ b/Src/Components/FBasicParticleContainer.hpp
@@ -103,7 +103,14 @@ public:
int getNbParticles() const{
return nbParticles;
}
-
+ /**
+ * @brief reset the number of particles
+ * @warning Only the number of particles is set to 0, the particles are still here.
+ */
+ void resetNumberOfParticles()
+ {
+ nbParticles = 0 ;
+ }
/**
* @brief getPositions
* @return a FReal*[3] to get access to the positions
@@ -116,7 +123,7 @@ public:
* @brief getWPositions
* @return get the position in write mode
*/
- FReal*const* getWPositions() {
+ FReal* const* getWPositions() {
return positions;
}
diff --git a/Src/Containers/FOctree.hpp b/Src/Containers/FOctree.hpp
index 733ca9525b1a55abf3f80fd630140db155f94bbf..bec09a4a97e4e2fff17f18b2dae04447572fd0a3 100755
--- a/Src/Containers/FOctree.hpp
+++ b/Src/Containers/FOctree.hpp
@@ -54,1123 +54,1183 @@
*/
template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass = FBasicBlockAllocator<CellClass> /*FListBlockAllocator<CellClass, 15>*/ >
class FOctree : public FNoCopyable {
- typedef FSubOctreeWithLeafs< CellClass , ContainerClass, LeafClass, CellAllocatorClass> SubOctreeWithLeaves;
- typedef FSubOctree< CellClass , ContainerClass, LeafClass, CellAllocatorClass> SubOctree;
-
- FAbstractSubOctree< CellClass , ContainerClass, LeafClass, CellAllocatorClass>* root; //< root suboctree
-
- FReal*const boxWidthAtLevel; //< to store the width of each boxs at all levels
-
- const int height; //< tree height
- const int subHeight; //< tree height
- const int leafIndex; //< index of leaf int array
-
- const FPoint boxCenter; //< the space system center
- const FPoint boxCorner; //< the space system corner (used to compute morton index)
-
- const FReal boxWidth; //< the space system width
-
-
- /**
- * Get morton index from a position for the leaf leavel
- * @param inPosition position to compute
- * @return the morton index
- */
- FTreeCoordinate getCoordinateFromPosition(const FPoint& inPosition) const {
- // box coordinate to host the particle
- FTreeCoordinate host;
- // position has to be relative to corner not center
- host.setX( getTreeCoordinate( inPosition.getX() - this->boxCorner.getX() ));
- host.setY( getTreeCoordinate( inPosition.getY() - this->boxCorner.getY() ));
- host.setZ( getTreeCoordinate( inPosition.getZ() - this->boxCorner.getZ() ));
- return host;
- }
-
- /**
- * Get the box number from a position
- * at a position POS with a leaf level box width of WI, the position is RELATIVE_TO_CORNER(POS)/WI
- * @param inRelativePosition a position from the corner of the box
- * @return the box num at the leaf level that contains inRelativePosition
- */
- int getTreeCoordinate(const FReal inRelativePosition) const {
- FAssertLF( (inRelativePosition >= 0 && inRelativePosition < this->boxWidth), "inRelativePosition : ",inRelativePosition );
- FAssertLF( inRelativePosition >= 0 && inRelativePosition < this->boxWidth, "Particle out of box" );
- const FReal indexFReal = inRelativePosition / this->boxWidthAtLevel[this->leafIndex];
- return static_cast<int>(indexFReal);
- }
+ typedef FOctree< CellClass , ContainerClass, LeafClass, CellAllocatorClass> OctreeType;
+ typedef FSubOctreeWithLeafs< CellClass , ContainerClass, LeafClass, CellAllocatorClass> SubOctreeWithLeaves;
+ typedef FSubOctree< CellClass , ContainerClass, LeafClass, CellAllocatorClass> SubOctree;
+
+ FAbstractSubOctree< CellClass , ContainerClass, LeafClass, CellAllocatorClass>* root; //< root suboctree
+
+ FReal*const boxWidthAtLevel; //< to store the width of each boxs at all levels
+
+ const int height; //< tree height
+ const int subHeight; //< tree height
+ const int leafIndex; //< index of leaf int array
+
+ const FPoint boxCenter; //< the space system center
+ const FPoint boxCorner; //< the space system corner (used to compute morton index)
+
+ const FReal boxWidth; //< the space system width
+
+
+ /**
+ * Get morton index from a position for the leaf leavel
+ * @param inPosition position to compute
+ * @return the morton index
+ */
+ FTreeCoordinate getCoordinateFromPosition(const FPoint& inPosition) const {
+ // box coordinate to host the particle
+ FTreeCoordinate host;
+ // position has to be relative to corner not center
+ host.setX( getTreeCoordinate( inPosition.getX() - this->boxCorner.getX() ));
+ host.setY( getTreeCoordinate( inPosition.getY() - this->boxCorner.getY() ));
+ host.setZ( getTreeCoordinate( inPosition.getZ() - this->boxCorner.getZ() ));
+ return host;
+ }
+
+ /**
+ * Get the box number from a position
+ * at a position POS with a leaf level box width of WI, the position is RELATIVE_TO_CORNER(POS)/WI
+ * @param inRelativePosition a position from the corner of the box
+ * @return the box num at the leaf level that contains inRelativePosition
+ */
+ int getTreeCoordinate(const FReal inRelativePosition) const {
+ FAssertLF( (inRelativePosition >= 0 && inRelativePosition < this->boxWidth), "inRelativePosition : ",inRelativePosition );
+ FAssertLF( inRelativePosition >= 0 && inRelativePosition < this->boxWidth, "Particle out of box" );
+ const FReal indexFReal = inRelativePosition / this->boxWidthAtLevel[this->leafIndex];
+ return static_cast<int>(indexFReal);
+ }
public:
- /**
- * Constructor
- * @param inHeight the octree height
- * @param inSubHeight the octree subheight
- * @param inBoxWidth box width for this simulation
- * @param inBoxCenter box center for this simulation
- */
- FOctree(const int inHeight, const int inSubHeight,
- const FReal inBoxWidth, const FPoint& inBoxCenter)
- : root(0), boxWidthAtLevel(new FReal[inHeight]),
- height(inHeight) , subHeight(inSubHeight), leafIndex(this->height-1),
- boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth)
- {
- FAssertLF(subHeight <= height - 1, "Subheight cannot be greater than height", __LINE__, __FILE__ );
- // Does we only need one suboctree?
- if(subHeight == height - 1){
- root = new FSubOctreeWithLeafs< CellClass , ContainerClass, LeafClass,CellAllocatorClass>(0, 0, this->subHeight, 1);
- }
- else {// if(subHeight < height - 1)
- root = new FSubOctree< CellClass , ContainerClass, LeafClass,CellAllocatorClass>(0, 0, this->subHeight, 1);
- }
-
- FReal tempWidth = this->boxWidth;
- // pre compute box width for each level
- for(int indexLevel = 0; indexLevel < this->height; ++indexLevel ){
- this->boxWidthAtLevel[indexLevel] = tempWidth;
- tempWidth /= FReal(2.0);
- }
- }
-
- /** Desctructor */
- virtual ~FOctree() {
- delete [] boxWidthAtLevel;
- delete root;
- }
-
- /** To get the tree height */
- int getHeight() const {
- return this->height;
- }
-
- /** To get the tree subheight */
- int getSubHeight() const{
- return this->subHeight;
- }
-
- /** To get the box width */
- FReal getBoxWidth() const{
- return this->boxWidth;
- }
-
- /** To get the center of the box */
- const FPoint& getBoxCenter() const{
- return this->boxCenter;
- }
-
- /** Count the number of cells per level,
- * it will iter on the tree to do that!
- */
- void getNbCellsPerLevel(int inNbCells[]){
- Iterator octreeIterator(this);
- octreeIterator.gotoBottomLeft();
-
- Iterator avoidGoLeft(octreeIterator);
-
- for(int idxLevel = height - 1 ; idxLevel > 0; --idxLevel ){
- int counter = 0;
- do{
- ++counter;
- } while(octreeIterator.moveRight());
- avoidGoLeft.moveUp();
- octreeIterator = avoidGoLeft;
- inNbCells[idxLevel] = counter;
- }
-
- inNbCells[0] = 0;
- }
-
- /**
- * Insert a particle on the tree
- * algorithm is :
- * Compute morton index for the particle
- * ask node to insert this particle
- * @param inParticle the particle to insert (must inherite from FAbstractParticle)
- */
- template<typename... Args>
- void insert(const FPoint& inParticlePosition, Args... args){
- const FTreeCoordinate host = getCoordinateFromPosition( inParticlePosition );
- const MortonIndex particleIndex = host.getMortonIndex(leafIndex);
- if(root->isLeafPart()){
- ((SubOctreeWithLeaves*)root)->insert( particleIndex, host, this->height, inParticlePosition, args... );
- }
- else{
- ((SubOctree*)root)->insert( particleIndex, host, this->height, inParticlePosition, args... );
- }
- }
-
- /** Remove a leaf from its morton index
- * @param indexToRemove the index of the leaf to remove
- */
- LeafClass* createLeaf(const MortonIndex indexToCreate ){
- const FTreeCoordinate host(indexToCreate,this->height-1);
- if(root->isLeafPart()){
- return ((SubOctreeWithLeaves*)root)->createLeaf( indexToCreate, host, this->height );
- }
- else{
- return ((SubOctree*)root)->createLeaf( indexToCreate, host, this->height );
- }
- }
-
- /** Remove a leaf from its morton index
- * @param indexToRemove the index of the leaf to remove
- */
- void removeLeaf(const MortonIndex indexToRemove ){
- root->removeLeaf( indexToRemove , this->height);
- }
-
- /**
- * Get a morton index from a real position
- * @param position a position to compute MI
- * @return the morton index
- */
- MortonIndex getMortonFromPosition(const FPoint& position) const {
- return getCoordinateFromPosition(position).getMortonIndex(leafIndex);
- }
-
- /**
- * The class works on suboctree. Most of the resources needed
- * are avaiblable by using FAbstractSubOctree. But when accessing
- * to the leaf we have to use FSubOctree or FSubOctreeWithLeafs
- * depending if we are working on the bottom of the tree.
- */
- union SubOctreeTypes {
- FAbstractSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree; //< Usual pointer to work
- FSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* middleTree; //< To access to sub-octree under
- FSubOctreeWithLeafs<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* leafTree;//< To access to particles lists
- };
-
- /**
- * This class is a const SubOctreeTypes
- */
- union SubOctreeTypesConst {
- const FAbstractSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree; //< Usual pointer to work
- const FSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* middleTree; //< To access to sub-octree under
- const FSubOctreeWithLeafs<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* leafTree;//< To access to particles lists
- };
-
- /**
- * This has to be used to iterate on an octree
- * It simply stores an pointer on a suboctree and moves to right/left/up/down.
- * Please refere to testOctreeIter file to see how it works.
- *
- * @code
- * FOctree<TestParticle, TestCell, NbLevels, NbSubLevels>::Iterator octreeIterator(&tree); <br>
- * octreeIterator.gotoBottomLeft(); <br>
- * for(int idx = 0 ; idx < NbLevels - 1; ++idx ){ <br>
- * do{ <br>
- * // ... <br>
- * } while(octreeIterator.moveRight()); <br>
- * octreeIterator.moveUp(); <br>
- * octreeIterator.gotoLeft(); <br>
- * } <br>
- * @endcode
- * Remark :
- * It uses the left right limit on each suboctree and their morton index.
- * Please have a look to the move functions to understand how the system is working.
- */
- class Iterator {
- SubOctreeTypes current; //< Current suboctree
-
- int currentLocalIndex; //< Current index (array position) in the current_suboctree.cells[ currentLocalLevel ]
- int currentLocalLevel; //< Current level in the current suboctree
-
- /**
- * To know what is the left limit on the current level on the current subtree
- * @return suboctree.left_limit >> 3 * diff(leaf_index, current_index).
- */
- static int TransposeIndex(const int indexInLeafLevel, const int distanceFromLeafLevel) {
- return indexInLeafLevel >> (3 * distanceFromLeafLevel);
- }
-
-
- public:
- /**
- * Constructor
- * @param inTarget the octree to iterate on
- * After building a iterator, this one is positioned at the level 0
- * of the root (level 1 of octree) at the left limit index
- */
- explicit Iterator(FOctree* const inTarget)
- : currentLocalIndex(0) , currentLocalLevel(0) {
- FAssertLF(inTarget, "Target for Octree::Iterator cannot be null", __LINE__, __FILE__);
- FAssertLF(inTarget->root->getRightLeafIndex() >= 0, "Octree seems to be empty, getRightLeafIndex == 0", __LINE__, __FILE__);
-
- // Start by the root
- this->current.tree = inTarget->root;
- // On the left limit
- this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1) );
- }
-
- Iterator() : currentLocalIndex(0),currentLocalLevel(0) {
- current.tree = 0;
- }
-
- /** Copy constructor
- * @param other source iterator to copy
- */
- Iterator(const Iterator& other){
- this->current = other.current ;
- this->currentLocalLevel = other.currentLocalLevel ;
- this->currentLocalIndex = other.currentLocalIndex ;
- }
-
- /** Copy operator
- * @param other source iterator to copy
- * @return this after copy
- */
- Iterator& operator=(const Iterator& other){
- this->current = other.current ;
- this->currentLocalLevel = other.currentLocalLevel ;
- this->currentLocalIndex = other.currentLocalIndex ;
- return *this;
- }
-
- /**
- * Move iterator to the top! (level 0 of root suboctree, level 1 of octree)
- * after this function : index = left limit at root level
- * the Algorithm is :
- * going to root suboctree
- * going to the first level and most left node
- */
- void gotoTop(){
- while(this->current.tree->hasParent()){
- this->current.tree = this->current.tree->getParent();
- }
- this->currentLocalLevel = 0;
- this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - 1) );
- }
-
- /**
- * Move iterator to the bottom left place
- * We are on a leaf a the most left node
- * the Algorithm is :
- * first go to top
- * then stay on the left and go downward
- */
- void gotoBottomLeft(){
- gotoTop();
- while(1) {
- this->currentLocalLevel = this->current.tree->getSubOctreeHeight() - 1;
- this->currentLocalIndex = this->current.tree->getLeftLeafIndex();
- if( isAtLeafLevel() ){
- return;
- }
- this->current.tree = this->current.middleTree->leafs( this->currentLocalIndex );
- }
- }
-
- /**
- * Move iterator to the left place at the same level
- * if needed we go on another suboctree but we stay on at the same level
- * the Algorithm is :
- * go to top
- * go downard until we are a the same level
- */
- void gotoLeft(){
- // Function variables
- const int currentLevel = level();
-
- // Goto root sutoctree
- while( this->current.tree->hasParent() ){
- this->current.tree = this->current.tree->getParent();
- }
-
- // Go down on the left until arriving on the same global level
- while( this->current.tree->getSubOctreeHeight() + this->current.tree->getSubOctreePosition() - 1 < currentLevel ) {
- this->current.tree = this->current.middleTree->leafs(this->current.tree->getLeftLeafIndex());
- }
-
- // Level still unchanged we only go to the left
- // The left limit on this tree at the level we want to stay
- this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1));
- }
-
- /**
- * Move iterator to the right place at the same level
- * if needed we go on another suboctree but we stay on at the same level
- * the Algorithm is :
- * go to top
- * go downard until we are a the same level
- */
- void gotoRight(){
- // Function variables
- const int currentLevel = level();
- // Goto root sutoctree
- while( this->current.tree->hasParent() ){
- this->current.tree = this->current.tree->getParent();
- }
- // Go down on the left until arriving on the same global level
- while( this->current.tree->getSubOctreeHeight() + this->current.tree->getSubOctreePosition() - 1 < currentLevel ) {
- this->current.tree = this->current.middleTree->leafs(this->current.tree->getRightLeafIndex());
- }
- // Level still unchanged we only go to the left
- // The left limit on this tree at the level we want to stay
- this->currentLocalIndex = TransposeIndex(this->current.tree->getRightLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1));
- }
-
- /**
- * Goto the next value on the right at the same level
- *
- * The algorithm here is :
- * As long as we are on the right limit, go to the parent suboctree
- * if we are on the root and on the right then return (there is no more data on the right)
- *
- * After that point we do not know where we are but we know that there is some data
- * on the right (without knowing our position!)
- *
- * We gotoNext on the brother to find an allocated cell (->)
- * for example if we are on index 2 we will look until 8 = 2 | 7 + 1
- * if we arrive a 8 without finding a cell we go upper and do the same
- * we know we will find something because we are not at the right limit
- *
- * We find an allocated cell.
- * We have to go down, we go on the left child of this cells
- * until : the current level if we did not have change the current suboctree
- * or : the leaf level
- *
- * In the second case, it meanse we need to change octree downward
- * but it is easy because we can use the left limit!
- *
- * @return true if we succeed to go to the right, else false
- */
- bool moveRight(){
- // Function variables
- SubOctreeTypes workingTree = this->current; // To cover parent other sub octre
- int workingLevel = this->currentLocalLevel; // To know where we are
- int workingIndex = this->currentLocalIndex; // To know where we are
-
- // -- First we have to go in a tree where we can move on the right --
- // Number of time we go into parent subtree
- int countUpward = 0;
- // We stop when we can right move or if there is no more parent (root)
- while( workingIndex == TransposeIndex(workingTree.tree->getRightLeafIndex(), (workingTree.tree->getSubOctreeHeight() - workingLevel - 1) )
- && workingTree.tree->hasParent() ){
- // Goto the leaf level into parent at current_tree.position_into_parent_array
- workingIndex = workingTree.tree->getIndexInParent();
- workingTree.tree = workingTree.tree->getParent();
- workingLevel = workingTree.tree->getSubOctreeHeight() - 1;
- // inc counter
- ++countUpward;
- }
-
- // Do we stop because we are on the root (means we cannot move right?)
- if( workingIndex < TransposeIndex(workingTree.tree->getRightLeafIndex(), (workingTree.tree->getSubOctreeHeight() - workingLevel - 1) ) ){
- // Move to the first right cell pointer(!)
- ++workingIndex;
-
- // Maybe it is null, but we know there is almost one cell on the right
- // we need to find it
- if( !workingTree.tree->cellsAt(workingLevel)[workingIndex] ){
- // While we are not on a allocated cell
- while( true ){
- // Test element on the right (test brothers)
- const int rightLimite = (workingIndex | 7) + 1;
- while( workingIndex < rightLimite && !workingTree.tree->cellsAt(workingLevel)[workingIndex]){
- ++workingIndex;
- }
- // Stop if we are on an allocated cell
- if( workingTree.tree->cellsAt(workingLevel)[workingIndex] ){
- break;
- }
- // Else go to the upper level
- --workingLevel;
- workingIndex >>= 3;
- }
- }
-
- // if wokring tree != current tree => working tree leafs level ; else current level
- const int objectiveLevel = (countUpward ? workingTree.tree->getSubOctreeHeight() - 1 : this->currentLocalLevel );
-
- // We need to go down as left as possible
- while( workingLevel != objectiveLevel ){
- ++workingLevel;
- workingIndex <<= 3;
- const int rightLimite = (workingIndex | 7); // not + 1 because if the 7th first are null it must be the 8th!
- while( workingIndex < rightLimite && !workingTree.tree->cellsAt(workingLevel)[workingIndex]){
- ++workingIndex;
- }
- }
-
- // Do we change from the current sub octree?
- if( countUpward ){
- // Then we simply need to go down the same number of time
- workingTree.tree = workingTree.middleTree->leafs(workingIndex);
- while( --countUpward ){
- workingTree.tree = workingTree.middleTree->leafs(workingTree.tree->getLeftLeafIndex());
- }
- // level did not change, simpli set octree and get left limit of this octree at the current level
- this->current = workingTree;
- this->currentLocalIndex = TransposeIndex(workingTree.tree->getLeftLeafIndex(), (workingTree.tree->getSubOctreeHeight() - this->currentLocalLevel - 1) );
- }
- else{
- // We are on the right tree
- this->currentLocalIndex = workingIndex;
- }
-
- return true;
- }
- return false;
- }
-
- /**
- * Move to the upper level
- * It may cause to change the suboctree we are working on
- * but we are using the same morton index >> 3
- * @return true if succeed
- */
- bool moveUp() {
- // It is on the top level?
- if( this->currentLocalLevel ){
- // No so simply go up
- --this->currentLocalLevel;
- this->currentLocalIndex >>= 3;
- }
- // Yes need to change suboctree
- else if( this->current.tree->hasParent() ){
- this->currentLocalIndex = this->current.tree->getIndexInParent();
- this->current.tree = this->current.tree->getParent();
- this->currentLocalLevel = this->current.tree->getSubOctreeHeight() - 1;
- }
- else{
- return false;
- }
- return true;
- }
-
- /**
- * Move down
- * It may cause to change the suboctree we are working on
- * We point on the first child found from left to right in the above
- * level
- * @return true if succeed
- */
- bool moveDown(){
- if( !isAtLeafLevel() ){
- // We are on the leaf of the current suboctree?
- if(this->currentLocalLevel + 1 == this->current.tree->getSubOctreeHeight()){
- // Yes change sub octree
- this->current.tree = this->current.middleTree->leafs(this->currentLocalIndex);
- this->currentLocalIndex = 0;
- this->currentLocalLevel = 0;
- }
- // No simply go down
- else{
- ++this->currentLocalLevel;
- this->currentLocalIndex <<= 3;
- }
- // Find the first allocated cell from left to right
- while(!this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]){
- ++this->currentLocalIndex;
- }
- return true;
- }
- return false;
- }
-
- /**
- * To know if we are not on the root level
- * @return true if we can move up
- */
- bool canProgressToUp() const {
- return this->currentLocalLevel || this->current.tree->hasParent();
- }
-
- /**
- * To know if we are not on the leafs level
- * @return true if we can move down
- */
- bool canProgressToDown() const {
- return !isAtLeafLevel();
- }
-
- /**
- * To know if we are on the leafs level
- * @return true if we are at the bottom of the tree
- */
- bool isAtLeafLevel() const {
- return this->current.tree->isLeafPart() && this->currentLocalLevel + 1 == this->current.tree->getSubOctreeHeight();
- }
-
- /**
- * To know the current level (not local but global)
- * @return the level in the entire octree
- */
- int level() const {
- return this->currentLocalLevel + this->current.tree->getSubOctreePosition();
- }
-
- /** Get the current pointed leaf
- * @return current leaf element
- */
- LeafClass* getCurrentLeaf() const {
- return this->current.leafTree->getLeaf(this->currentLocalIndex);
- }
-
- /** To access the current particles list
- * You have to be at the leaf level to call this function!
- * @return current element list
- */
- ContainerClass* getCurrentListSrc() const {
- return this->current.leafTree->getLeafSrc(this->currentLocalIndex);
- }
-
- /** To access the current particles list
- * You have to be at the leaf level to call this function!
- * @return current element list
- */
- ContainerClass* getCurrentListTargets() const {
- return this->current.leafTree->getLeafTargets(this->currentLocalIndex);
- }
-
- /** Get the current pointed cell
- * @return current cell element
- */
- CellClass* getCurrentCell() const {
- return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex];
- }
-
- /** Get the child of the current cell
- * This function return an array of CellClass (array size = 8)
- * User has to test each case to know if there is a cell
- * @return the child array
- */
- CellClass** getCurrentChild() const {
- // are we at the bottom of the suboctree
- if(this->current.tree->getSubOctreeHeight() - 1 == this->currentLocalLevel ){
- // then return first level of the suboctree under
- return &this->current.middleTree->leafs(this->currentLocalIndex)->cellsAt(0)[0];
- }
- else{
- // else simply return the array at the right position
- return &this->current.tree->cellsAt(this->currentLocalLevel + 1)[this->currentLocalIndex << 3];
- }
- }
-
- /** Get the part of array that contains all the pointers
- *
- */
- CellClass** getCurrentBox() const {
- return &this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex & ~7];
- }
-
- /** Get the Morton index of the current cell pointed by the iterator
- * @return The global morton index
- * <code>iter.getCurrentGlobalIndex();<br>
- * // is equivalent to :<br>
- * iter.getCurrentCell()->getMortonIndex();</code>
- */
- MortonIndex getCurrentGlobalIndex() const{
- return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]->getMortonIndex();
- }
-
- /** To get the tree coordinate of the current working cell
- *
- */
- const FTreeCoordinate& getCurrentGlobalCoordinate() const{
- return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]->getCoordinate();
- }
-
- };
-
- // To be able to access octree root & data
- 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 find a missing cell or the right cell
- */
- CellClass* getCell(const MortonIndex inIndex, const int inLevel) const{
- SubOctreeTypesConst workingTree;
- workingTree.tree = &this->root;
- const MortonIndex treeSubLeafMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
-
- // Find the suboctree a the correct level
- while(inLevel >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
- // compute the leaf index
- const MortonIndex fullIndex = inIndex >> (3 * (inLevel + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) ) );
- // point to next suboctree
- workingTree.tree = workingTree.middleTree->leafs(treeSubLeafMask & fullIndex);
- if(!workingTree.tree) return 0;
- }
-
- // compute correct index in the array
- const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (inLevel + 1 - workingTree.tree->getSubOctreePosition()) ));
- return workingTree.tree->cellsAt(inLevel - workingTree.tree->getSubOctreePosition())[treeLeafMask & inIndex];
- }
-
-
- /** This function fill a array with the neighbors of a cell
- * it does not put the brothers in the array (brothers are cells
- * at the same level with the same parent) because they are of course
- * direct neighbors.
- * There is a maximum of 26 (3*3*3-1) direct neighbors
- * // Take the neighbors != brothers
- * CellClass* directNeighbors[26];
- * const int nbDirectNeighbors = getNeighborsNoBrothers(directNeighbors,inIndex,inLevel);
- * @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 getNeighborsNoBrothers(CellClass* inNeighbors[26], const MortonIndex inIndex, const int inLevel) const {
- FTreeCoordinate center;
- center.setPositionFromMorton(inIndex, inLevel);
-
- const int boxLimite = FMath::pow2(inLevel);
-
- int idxNeighbors = 0;
-
- // We test all cells around
- for(int idxX = -1 ; idxX <= 1 ; ++idxX){
- if(!FMath::Between(center.getX() + idxX,0,boxLimite)) continue;
-
- for(int idxY = -1 ; idxY <= 1 ; ++idxY){
- if(!FMath::Between(center.getY() + idxY,0,boxLimite)) continue;
-
- for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
- if(!FMath::Between(center.getZ() + idxZ,0,boxLimite)) continue;
-
- // 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(inLevel);
- // if not a brother
- if( mortonOther>>3 != inIndex>>3 ){
- // get cell
- CellClass* const cell = getCell(mortonOther, inLevel);
- // add to list if not null
- if(cell) inNeighbors[idxNeighbors++] = cell;
- }
- }
- }
- }
- }
-
- return idxNeighbors;
- }
-
-
- /** This function return an adresse of cell array from a morton index and a level
- *
- * @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)
- *
- */
- CellClass** getCellPt(const MortonIndex inIndex, const int inLevel) const{
- SubOctreeTypesConst workingTree;
- workingTree.tree = this->root;
-
- const MortonIndex treeMiddleMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
-
- // Find the suboctree a the correct level
- while(inLevel >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
- // compute the leaf index
- const MortonIndex fullIndex = inIndex >> 3 * (inLevel + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) );
- // point to next suboctree
- workingTree.tree = workingTree.middleTree->leafs(int(treeMiddleMask & fullIndex));
- if(!workingTree.tree) return 0;
- }
-
- // Be sure there is a parent allocated
- const int levelInTree = inLevel - workingTree.tree->getSubOctreePosition();
- if( levelInTree && !workingTree.tree->cellsAt(levelInTree - 1)[~(~0x00LL << (3 * levelInTree )) & (inIndex>>3)]){
- return 0;
- }
-
- // compute correct index in the array and return the @ in array
- const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (levelInTree + 1 ) ));
- return &workingTree.tree->cellsAt(levelInTree)[treeLeafMask & inIndex];
- }
-
-
- /** This function fill an array with the distant neighbors of a cell
- * @param inNeighbors the array to store the elements
- * @param inNeighborsIndex the array to store morton index of the neighbors
- * @param inIndex the index of the element we want the neighbors
- * @param inLevel the level of the element
- * @return the number of neighbors
- */
- int getInteractionNeighbors(const CellClass* inNeighbors[343],
- const FTreeCoordinate& workingCell,
- const int inLevel) const{
- // reset
- memset(inNeighbors, 0, sizeof(CellClass*) * 343);
-
- // Then take each child of the parent's neighbors if not in directNeighbors
- // Father coordinate
- const FTreeCoordinate parentCell(workingCell.getX()>>1,workingCell.getY()>>1,workingCell.getZ()>>1);
-
- // Limite at parent level number of box (split by 2 by level)
- const int boxLimite = FMath::pow2(inLevel-1);
-
- int idxNeighbors = 0;
- // We test all cells around
- for(int idxX = -1 ; idxX <= 1 ; ++idxX){
- if(!FMath::Between(parentCell.getX() + idxX,0,boxLimite)) continue;
-
- for(int idxY = -1 ; idxY <= 1 ; ++idxY){
- if(!FMath::Between(parentCell.getY() + idxY,0,boxLimite)) continue;
-
- for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
- if(!FMath::Between(parentCell.getZ() + idxZ,0,boxLimite)) continue;
-
- // if we are not on the current cell
- if( idxX || idxY || idxZ ){
- const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
- const MortonIndex mortonOtherParent = otherParent.getMortonIndex(inLevel-1) << 3;
- // Get child
- CellClass** const cells = getCellPt(mortonOtherParent, inLevel);
-
- // If there is one or more child
- if(cells){
- // For each child
- for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
- if(cells[idxCousin]){
- const int xdiff = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - workingCell.getX();
- const int ydiff = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - workingCell.getY();
- const int zdiff = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
-
- // Test if it is a direct neighbor
- if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
- // add to neighbors
- inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
- ++idxNeighbors;
- }
- }
- }
- }
- }
- }
- }
- }
-
- return idxNeighbors;
- }
-
-
- /** This function fill an array with the distant neighbors of a cell
- * it respects the periodic condition and will give the relative distance
- * between the working cell and the neighbors
- * @param inNeighbors the array to store the elements
- * @param inRelativePosition the array to store the relative position of the neighbors
- * @param workingCell the index of the element we want the neighbors
- * @param inLevel the level of the element
- * @return the number of neighbors
- */
- int getPeriodicInteractionNeighbors(const CellClass* inNeighbors[343],
- const FTreeCoordinate& workingCell,
- const int inLevel, const int inDirection) const{
-
- // Then take each child of the parent's neighbors if not in directNeighbors
- // Father coordinate
- const FTreeCoordinate parentCell(workingCell.getX()>>1,workingCell.getY()>>1,workingCell.getZ()>>1);
-
- // Limite at parent level number of box (split by 2 by level)
- const int boxLimite = FMath::pow2(inLevel-1);
-
- // This is not on a border we can use normal interaction list method
- if( !(parentCell.getX() == 0 || parentCell.getY() == 0 || parentCell.getZ() == 0 ||
- parentCell.getX() == boxLimite - 1 || parentCell.getY() == boxLimite - 1 || parentCell.getZ() == boxLimite - 1 ) ) {
- return getInteractionNeighbors( inNeighbors, workingCell, inLevel);
- }
- else{
- // reset
- memset(inNeighbors, 0, sizeof(CellClass*) * 343);
-
- const int startX = (TestPeriodicCondition(inDirection, DirMinusX) || parentCell.getX() != 0 ?-1:0);
- const int endX = (TestPeriodicCondition(inDirection, DirPlusX) || parentCell.getX() != boxLimite - 1 ?1:0);
- const int startY = (TestPeriodicCondition(inDirection, DirMinusY) || parentCell.getY() != 0 ?-1:0);
- const int endY = (TestPeriodicCondition(inDirection, DirPlusY) || parentCell.getY() != boxLimite - 1 ?1:0);
- const int startZ = (TestPeriodicCondition(inDirection, DirMinusZ) || parentCell.getZ() != 0 ?-1:0);
- const int endZ = (TestPeriodicCondition(inDirection, DirPlusZ) || parentCell.getZ() != boxLimite - 1 ?1:0);
-
- int idxNeighbors = 0;
- // We test all cells around
- for(int idxX = startX ; idxX <= endX ; ++idxX){
- for(int idxY = startY ; idxY <= endY ; ++idxY){
- for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
- // if we are not on the current cell
- if( idxX || idxY || idxZ ){
-
- const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
- FTreeCoordinate otherParentInBox(otherParent);
- // periodic
- if( otherParentInBox.getX() < 0 ){
- otherParentInBox.setX( otherParentInBox.getX() + boxLimite );
- }
- else if( boxLimite <= otherParentInBox.getX() ){
- otherParentInBox.setX( otherParentInBox.getX() - boxLimite );
- }
-
- if( otherParentInBox.getY() < 0 ){
- otherParentInBox.setY( otherParentInBox.getY() + boxLimite );
- }
- else if( boxLimite <= otherParentInBox.getY() ){
- otherParentInBox.setY( otherParentInBox.getY() - boxLimite );
- }
-
- if( otherParentInBox.getZ() < 0 ){
- otherParentInBox.setZ( otherParentInBox.getZ() + boxLimite );
- }
- else if( boxLimite <= otherParentInBox.getZ() ){
- otherParentInBox.setZ( otherParentInBox.getZ() - boxLimite );
- }
-
-
- const MortonIndex mortonOtherParent = otherParentInBox.getMortonIndex(inLevel-1) << 3;
- // Get child
- CellClass** const cells = getCellPt(mortonOtherParent, inLevel);
-
- // If there is one or more child
- if(cells){
- // For each child
- for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
- if(cells[idxCousin]){
- const int xdiff = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - workingCell.getX();
- const int ydiff = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - workingCell.getY();
- const int zdiff = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
-
- // Test if it is a direct neighbor
- if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
- // add to neighbors
- inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
- ++idxNeighbors;
- }
- }
- }
- }
- }
- }
- }
- }
-
- return idxNeighbors;
- }
- }
-
-
- /** 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)
- *
- */
- ContainerClass* getLeafSrc(const MortonIndex inIndex){
- SubOctreeTypes workingTree;
- workingTree.tree = this->root;
- const MortonIndex treeSubLeafMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
-
- // Find the suboctree a the correct level
- while(leafIndex >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
- // compute the leaf index
- const MortonIndex fullIndex = inIndex >> (3 * (leafIndex + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) ) );
- // point to next suboctree
- workingTree.tree = workingTree.middleTree->leafs(int(treeSubLeafMask & fullIndex));
- if(!workingTree.tree) return 0;
- }
-
- // compute correct index in the array
- const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (leafIndex + 1 - workingTree.tree->getSubOctreePosition()) ));
- return workingTree.leafTree->getLeafSrc(int(treeLeafMask & inIndex));
- }
-
- /** 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(ContainerClass* inNeighbors[27], const FTreeCoordinate& center, const int inLevel){
- memset( inNeighbors, 0 , 27 * sizeof(ContainerClass*));
- const int boxLimite = FMath::pow2(inLevel);
-
- int idxNeighbors = 0;
-
- // We test all cells around
- for(int idxX = -1 ; idxX <= 1 ; ++idxX){
- if(!FMath::Between(center.getX() + idxX,0,boxLimite)) continue;
-
- for(int idxY = -1 ; idxY <= 1 ; ++idxY){
- if(!FMath::Between(center.getY() + idxY,0,boxLimite)) continue;
-
- for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
- if(!FMath::Between(center.getZ() + idxZ,0,boxLimite)) continue;
-
- // 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(inLevel);
- // get cell
- ContainerClass* const leaf = getLeafSrc(mortonOther);
- // add to list if not null
- if(leaf){
- inNeighbors[(((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1] = leaf;
- ++idxNeighbors;
- }
- }
- }
- }
- }
-
- return idxNeighbors;
- }
-
-
- /** 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 getPeriodicLeafsNeighbors(ContainerClass* inNeighbors[27], FTreeCoordinate inOffsets[27], bool*const isPeriodic,
- const FTreeCoordinate& center, const int inLevel, const int inDirection){
-
- const int boxLimite = FMath::pow2(inLevel);
-
- if( center.getX() != 0 && center.getY() != 0 && center.getZ() != 0 &&
- center.getX() != boxLimite - 1 && center.getY() != boxLimite - 1 && center.getZ() != boxLimite - 1 ){
- (*isPeriodic) = false;
- return getLeafsNeighbors(inNeighbors, center, inLevel);
- }
-
- (*isPeriodic) = true;
- memset(inNeighbors , 0 , 27 * sizeof(ContainerClass*));
- int idxNeighbors = 0;
-
- const int startX = (TestPeriodicCondition(inDirection, DirMinusX) || center.getX() != 0 ?-1:0);
- const int endX = (TestPeriodicCondition(inDirection, DirPlusX) || center.getX() != boxLimite - 1 ?1:0);
- const int startY = (TestPeriodicCondition(inDirection, DirMinusY) || center.getY() != 0 ?-1:0);
- const int endY = (TestPeriodicCondition(inDirection, DirPlusY) || center.getY() != boxLimite - 1 ?1:0);
- const int startZ = (TestPeriodicCondition(inDirection, DirMinusZ) || center.getZ() != 0 ?-1:0);
- const int endZ = (TestPeriodicCondition(inDirection, DirPlusZ) || center.getZ() != boxLimite - 1 ?1:0);
-
- // We test all cells around
- for(int idxX = startX ; idxX <= endX ; ++idxX){
- for(int idxY = startY ; idxY <= endY ; ++idxY){
- for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
- // if we are not on the current cell
- if( idxX || idxY || idxZ ){
- FTreeCoordinate other(center.getX() + idxX,center.getY() + idxY,center.getZ() + idxZ);
- int xoffset = 0, yoffset = 0, zoffset = 0;
-
- if( other.getX() < 0 ){
- other.setX( other.getX() + boxLimite );
- xoffset = -1;
- }
- else if( boxLimite <= other.getX() ){
- other.setX( other.getX() - boxLimite );
- xoffset = 1;
- }
- if( other.getY() < 0 ){
- other.setY( other.getY() + boxLimite );
- yoffset = -1;
- }
- else if( boxLimite <= other.getY() ){
- other.setY( other.getY() - boxLimite );
- yoffset = 1;
- }
- if( other.getZ() < 0 ){
- other.setZ( other.getZ() + boxLimite );
- zoffset = -1;
- }
- else if( boxLimite <= other.getZ() ){
- other.setZ( other.getZ() - boxLimite );
- zoffset = 1;
- }
-
- const MortonIndex mortonOther = other.getMortonIndex(inLevel);
- // get cell
- ContainerClass* const leaf = getLeafSrc(mortonOther);
- // add to list if not null
- if(leaf){
- const int index = (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1;
- inNeighbors[index] = leaf;
- inOffsets[index].setPosition(xoffset,yoffset,zoffset);
- ++idxNeighbors;
- }
- }
- }
- }
- }
-
- return idxNeighbors;
- }
-
-
-
- /////////////////////////////////////////////////////////
- // Lambda function to apply to all member
- /////////////////////////////////////////////////////////
-
- /**
- * @brief forEachLeaf iterate on the leaf and apply the function
- * @param function
- */
- void forEachLeaf(std::function<void(LeafClass*)> function){
- Iterator octreeIterator(this);
- octreeIterator.gotoBottomLeft();
-
- do{
- function(octreeIterator.getCurrentLeaf());
- } while(octreeIterator.moveRight());
- }
-
- /**
- * @brief forEachLeaf iterate on the cell and apply the function
- * @param function
- */
- void forEachCell(std::function<void(CellClass*)> function){
- Iterator octreeIterator(this);
- octreeIterator.gotoBottomLeft();
-
- Iterator avoidGoLeft(octreeIterator);
-
- for(int idx = this->height-1 ; idx >= 1 ; --idx ){
- do{
- function(octreeIterator.getCurrentCell());
- } while(octreeIterator.moveRight());
- avoidGoLeft.moveUp();
- octreeIterator = avoidGoLeft;
- }
- }
-
- /**
- * @brief forEachLeaf iterate on the cell and apply the function
- * @param function
- */
- void forEachCellWithLevel(std::function<void(CellClass*,const int)> function){
- Iterator octreeIterator(this);
- octreeIterator.gotoBottomLeft();
-
- Iterator avoidGoLeft(octreeIterator);
-
- for(int idx = this->height-1 ; idx >= 1 ; --idx ){
- do{
- function(octreeIterator.getCurrentCell(),idx);
- } while(octreeIterator.moveRight());
- avoidGoLeft.moveUp();
- octreeIterator = avoidGoLeft;
- }
- }
-
- /**
- * @brief forEachLeaf iterate on the cell and apply the function
- * @param function
- */
- void forEachCellLeaf(std::function<void(CellClass*,LeafClass*)> function){
- Iterator octreeIterator(this);
- octreeIterator.gotoBottomLeft();
-
- do{
- function(octreeIterator.getCurrentCell(),octreeIterator.getCurrentLeaf());
- } while(octreeIterator.moveRight());
- }
+ /**
+ * Constructor
+ * @param inHeight the octree height
+ * @param inSubHeight the octree subheight
+ * @param inBoxWidth box width for this simulation
+ * @param inBoxCenter box center for this simulation
+ */
+ FOctree(const int inHeight, const int inSubHeight,
+ const FReal inBoxWidth, const FPoint& inBoxCenter)
+: root(0), boxWidthAtLevel(new FReal[inHeight]),
+ height(inHeight) , subHeight(inSubHeight), leafIndex(this->height-1),
+ boxCenter(inBoxCenter), boxCorner(inBoxCenter,-(inBoxWidth/2)), boxWidth(inBoxWidth)
+{
+ FAssertLF(subHeight <= height - 1, "Subheight cannot be greater than height", __LINE__, __FILE__ );
+ // Does we only need one suboctree?
+ if(subHeight == height - 1){
+ root = new FSubOctreeWithLeafs< CellClass , ContainerClass, LeafClass,CellAllocatorClass>(0, 0, this->subHeight, 1);
+ }
+ else {// if(subHeight < height - 1)
+ root = new FSubOctree< CellClass , ContainerClass, LeafClass,CellAllocatorClass>(0, 0, this->subHeight, 1);
+ }
+
+ FReal tempWidth = this->boxWidth;
+ // pre compute box width for each level
+ for(int indexLevel = 0; indexLevel < this->height; ++indexLevel ){
+ this->boxWidthAtLevel[indexLevel] = tempWidth;
+ tempWidth /= FReal(2.0);
+ }
+}
+
+ /** Desctructor */
+ virtual ~FOctree() {
+ delete [] boxWidthAtLevel;
+ delete root;
+ }
+
+ /** To get the tree height */
+ int getHeight() const {
+ return this->height;
+ }
+
+ /** To get the tree subheight */
+ int getSubHeight() const{
+ return this->subHeight;
+ }
+
+ /** To get the box width */
+ FReal getBoxWidth() const{
+ return this->boxWidth;
+ }
+
+ /** To get the center of the box */
+ const FPoint& getBoxCenter() const{
+ return this->boxCenter;
+ }
+
+ /** Count the number of cells per level,
+ * it will iter on the tree to do that!
+ */
+ void getNbCellsPerLevel(int inNbCells[]){
+ Iterator octreeIterator(this);
+ octreeIterator.gotoBottomLeft();
+
+ Iterator avoidGoLeft(octreeIterator);
+
+ for(int idxLevel = height - 1 ; idxLevel > 0; --idxLevel ){
+ int counter = 0;
+ do{
+ ++counter;
+ } while(octreeIterator.moveRight());
+ avoidGoLeft.moveUp();
+ octreeIterator = avoidGoLeft;
+ inNbCells[idxLevel] = counter;
+ }
+
+ inNbCells[0] = 0;
+ }
+
+ /**
+ * Insert a particle on the tree
+ * algorithm is :
+ * Compute morton index for the particle
+ * ask node to insert this particle
+ * @param inParticle the particle to insert (must inherite from FAbstractParticle)
+ */
+ template<typename... Args>
+ void insert(const FPoint& inParticlePosition, Args... args){
+ const FTreeCoordinate host = getCoordinateFromPosition( inParticlePosition );
+ const MortonIndex particleIndex = host.getMortonIndex(leafIndex);
+ if(root->isLeafPart()){
+ ((SubOctreeWithLeaves*)root)->insert( particleIndex, host, this->height, inParticlePosition, args... );
+ }
+ else{
+ ((SubOctree*)root)->insert( particleIndex, host, this->height, inParticlePosition, args... );
+ }
+ }
+
+ /** Remove a leaf from its morton index
+ * @param indexToRemove the index of the leaf to remove
+ */
+ LeafClass* createLeaf(const MortonIndex indexToCreate ){
+ const FTreeCoordinate host(indexToCreate,this->height-1);
+ if(root->isLeafPart()){
+ return ((SubOctreeWithLeaves*)root)->createLeaf( indexToCreate, host, this->height );
+ }
+ else{
+ return ((SubOctree*)root)->createLeaf( indexToCreate, host, this->height );
+ }
+ }
+
+ /** Remove a leaf from its morton index
+ * @param indexToRemove the index of the leaf to remove
+ */
+ void removeLeaf(const MortonIndex indexToRemove ){
+ root->removeLeaf( indexToRemove , this->height);
+ }
+
+ /**
+ * Get a morton index from a real position
+ * @param position a position to compute MI
+ * @return the morton index
+ */
+ MortonIndex getMortonFromPosition(const FPoint& position) const {
+ return getCoordinateFromPosition(position).getMortonIndex(leafIndex);
+ }
+
+ /**
+ * The class works on suboctree. Most of the resources needed
+ * are avaiblable by using FAbstractSubOctree. But when accessing
+ * to the leaf we have to use FSubOctree or FSubOctreeWithLeafs
+ * depending if we are working on the bottom of the tree.
+ */
+ union SubOctreeTypes {
+ FAbstractSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree; //< Usual pointer to work
+ FSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* middleTree; //< To access to sub-octree under
+ FSubOctreeWithLeafs<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* leafTree;//< To access to particles lists
+ };
+
+ /**
+ * This class is a const SubOctreeTypes
+ */
+ union SubOctreeTypesConst {
+ const FAbstractSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree; //< Usual pointer to work
+ const FSubOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* middleTree; //< To access to sub-octree under
+ const FSubOctreeWithLeafs<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* leafTree;//< To access to particles lists
+ };
+
+ /**
+ * This has to be used to iterate on an octree
+ * It simply stores an pointer on a suboctree and moves to right/left/up/down.
+ * Please refere to testOctreeIter file to see how it works.
+ *
+ * @code
+ * FOctree<TestParticle, TestCell, NbLevels, NbSubLevels>::Iterator octreeIterator(&tree); <br>
+ * octreeIterator.gotoBottomLeft(); <br>
+ * for(int idx = 0 ; idx < NbLevels - 1; ++idx ){ <br>
+ * do{ <br>
+ * // ... <br>
+ * } while(octreeIterator.moveRight()); <br>
+ * octreeIterator.moveUp(); <br>
+ * octreeIterator.gotoLeft(); <br>
+ * } <br>
+ * @endcode
+ * Remark :
+ * It uses the left right limit on each suboctree and their morton index.
+ * Please have a look to the move functions to understand how the system is working.
+ */
+ class Iterator {
+ SubOctreeTypes current; //< Current suboctree
+
+ int currentLocalIndex; //< Current index (array position) in the current_suboctree.cells[ currentLocalLevel ]
+ int currentLocalLevel; //< Current level in the current suboctree
+
+ /**
+ * To know what is the left limit on the current level on the current subtree
+ * @return suboctree.left_limit >> 3 * diff(leaf_index, current_index).
+ */
+ static int TransposeIndex(const int indexInLeafLevel, const int distanceFromLeafLevel) {
+ return indexInLeafLevel >> (3 * distanceFromLeafLevel);
+ }
+
+
+ public:
+ /**
+ * Constructor
+ * @param inTarget the octree to iterate on
+ * After building a iterator, this one is positioned at the level 0
+ * of the root (level 1 of octree) at the left limit index
+ */
+ explicit Iterator(FOctree* const inTarget)
+ : currentLocalIndex(0) , currentLocalLevel(0) {
+ FAssertLF(inTarget, "Target for Octree::Iterator cannot be null", __LINE__, __FILE__);
+ FAssertLF(inTarget->root->getRightLeafIndex() >= 0, "Octree seems to be empty, getRightLeafIndex == 0", __LINE__, __FILE__);
+
+ // Start by the root
+ this->current.tree = inTarget->root;
+ // On the left limit
+ this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1) );
+ }
+
+ Iterator() : currentLocalIndex(0),currentLocalLevel(0) {
+ current.tree = 0;
+ }
+
+ /** Copy constructor
+ * @param other source iterator to copy
+ */
+ Iterator(const Iterator& other){
+ this->current = other.current ;
+ this->currentLocalLevel = other.currentLocalLevel ;
+ this->currentLocalIndex = other.currentLocalIndex ;
+ }
+
+ /** Copy operator
+ * @param other source iterator to copy
+ * @return this after copy
+ */
+ Iterator& operator=(const Iterator& other){
+ this->current = other.current ;
+ this->currentLocalLevel = other.currentLocalLevel ;
+ this->currentLocalIndex = other.currentLocalIndex ;
+ return *this;
+ }
+
+ /**
+ * Move iterator to the top! (level 0 of root suboctree, level 1 of octree)
+ * after this function : index = left limit at root level
+ * the Algorithm is :
+ * going to root suboctree
+ * going to the first level and most left node
+ */
+ void gotoTop(){
+ while(this->current.tree->hasParent()){
+ this->current.tree = this->current.tree->getParent();
+ }
+ this->currentLocalLevel = 0;
+ this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - 1) );
+ }
+
+ /**
+ * Move iterator to the bottom left place
+ * We are on a leaf a the most left node
+ * the Algorithm is :
+ * first go to top
+ * then stay on the left and go downward
+ */
+ void gotoBottomLeft(){
+ gotoTop();
+ while(1) {
+ this->currentLocalLevel = this->current.tree->getSubOctreeHeight() - 1;
+ this->currentLocalIndex = this->current.tree->getLeftLeafIndex();
+ if( isAtLeafLevel() ){
+ return;
+ }
+ this->current.tree = this->current.middleTree->leafs( this->currentLocalIndex );
+ }
+ }
+
+ /**
+ * Move iterator to the left place at the same level
+ * if needed we go on another suboctree but we stay on at the same level
+ * the Algorithm is :
+ * go to top
+ * go downard until we are a the same level
+ */
+ void gotoLeft(){
+ // Function variables
+ const int currentLevel = level();
+
+ // Goto root sutoctree
+ while( this->current.tree->hasParent() ){
+ this->current.tree = this->current.tree->getParent();
+ }
+
+ // Go down on the left until arriving on the same global level
+ while( this->current.tree->getSubOctreeHeight() + this->current.tree->getSubOctreePosition() - 1 < currentLevel ) {
+ this->current.tree = this->current.middleTree->leafs(this->current.tree->getLeftLeafIndex());
+ }
+
+ // Level still unchanged we only go to the left
+ // The left limit on this tree at the level we want to stay
+ this->currentLocalIndex = TransposeIndex(this->current.tree->getLeftLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1));
+ }
+
+ /**
+ * Move iterator to the right place at the same level
+ * if needed we go on another suboctree but we stay on at the same level
+ * the Algorithm is :
+ * go to top
+ * go downard until we are a the same level
+ */
+ void gotoRight(){
+ // Function variables
+ const int currentLevel = level();
+ // Goto root sutoctree
+ while( this->current.tree->hasParent() ){
+ this->current.tree = this->current.tree->getParent();
+ }
+ // Go down on the left until arriving on the same global level
+ while( this->current.tree->getSubOctreeHeight() + this->current.tree->getSubOctreePosition() - 1 < currentLevel ) {
+ this->current.tree = this->current.middleTree->leafs(this->current.tree->getRightLeafIndex());
+ }
+ // Level still unchanged we only go to the left
+ // The left limit on this tree at the level we want to stay
+ this->currentLocalIndex = TransposeIndex(this->current.tree->getRightLeafIndex(), (this->current.tree->getSubOctreeHeight() - this->currentLocalLevel - 1));
+ }
+
+ /**
+ * Goto the next value on the right at the same level
+ *
+ * The algorithm here is :
+ * As long as we are on the right limit, go to the parent suboctree
+ * if we are on the root and on the right then return (there is no more data on the right)
+ *
+ * After that point we do not know where we are but we know that there is some data
+ * on the right (without knowing our position!)
+ *
+ * We gotoNext on the brother to find an allocated cell (->)
+ * for example if we are on index 2 we will look until 8 = 2 | 7 + 1
+ * if we arrive a 8 without finding a cell we go upper and do the same
+ * we know we will find something because we are not at the right limit
+ *
+ * We find an allocated cell.
+ * We have to go down, we go on the left child of this cells
+ * until : the current level if we did not have change the current suboctree
+ * or : the leaf level
+ *
+ * In the second case, it meanse we need to change octree downward
+ * but it is easy because we can use the left limit!
+ *
+ * @return true if we succeed to go to the right, else false
+ */
+ bool moveRight(){
+ // Function variables
+ SubOctreeTypes workingTree = this->current; // To cover parent other sub octre
+ int workingLevel = this->currentLocalLevel; // To know where we are
+ int workingIndex = this->currentLocalIndex; // To know where we are
+
+ // -- First we have to go in a tree where we can move on the right --
+ // Number of time we go into parent subtree
+ int countUpward = 0;
+ // We stop when we can right move or if there is no more parent (root)
+ while( workingIndex == TransposeIndex(workingTree.tree->getRightLeafIndex(), (workingTree.tree->getSubOctreeHeight() - workingLevel - 1) )
+ && workingTree.tree->hasParent() ){
+ // Goto the leaf level into parent at current_tree.position_into_parent_array
+ workingIndex = workingTree.tree->getIndexInParent();
+ workingTree.tree = workingTree.tree->getParent();
+ workingLevel = workingTree.tree->getSubOctreeHeight() - 1;
+ // inc counter
+ ++countUpward;
+ }
+
+ // Do we stop because we are on the root (means we cannot move right?)
+ if( workingIndex < TransposeIndex(workingTree.tree->getRightLeafIndex(), (workingTree.tree->getSubOctreeHeight() - workingLevel - 1) ) ){
+ // Move to the first right cell pointer(!)
+ ++workingIndex;
+
+ // Maybe it is null, but we know there is almost one cell on the right
+ // we need to find it
+ if( !workingTree.tree->cellsAt(workingLevel)[workingIndex] ){
+ // While we are not on a allocated cell
+ while( true ){
+ // Test element on the right (test brothers)
+ const int rightLimite = (workingIndex | 7) + 1;
+ while( workingIndex < rightLimite && !workingTree.tree->cellsAt(workingLevel)[workingIndex]){
+ ++workingIndex;
+ }
+ // Stop if we are on an allocated cell
+ if( workingTree.tree->cellsAt(workingLevel)[workingIndex] ){
+ break;
+ }
+ // Else go to the upper level
+ --workingLevel;
+ workingIndex >>= 3;
+ }
+ }
+
+ // if wokring tree != current tree => working tree leafs level ; else current level
+ const int objectiveLevel = (countUpward ? workingTree.tree->getSubOctreeHeight() - 1 : this->currentLocalLevel );
+
+ // We need to go down as left as possible
+ while( workingLevel != objectiveLevel ){
+ ++workingLevel;
+ workingIndex <<= 3;
+ const int rightLimite = (workingIndex | 7); // not + 1 because if the 7th first are null it must be the 8th!
+ while( workingIndex < rightLimite && !workingTree.tree->cellsAt(workingLevel)[workingIndex]){
+ ++workingIndex;
+ }
+ }
+
+ // Do we change from the current sub octree?
+ if( countUpward ){
+ // Then we simply need to go down the same number of time
+ workingTree.tree = workingTree.middleTree->leafs(workingIndex);
+ while( --countUpward ){
+ workingTree.tree = workingTree.middleTree->leafs(workingTree.tree->getLeftLeafIndex());
+ }
+ // level did not change, simpli set octree and get left limit of this octree at the current level
+ this->current = workingTree;
+ this->currentLocalIndex = TransposeIndex(workingTree.tree->getLeftLeafIndex(), (workingTree.tree->getSubOctreeHeight() - this->currentLocalLevel - 1) );
+ }
+ else{
+ // We are on the right tree
+ this->currentLocalIndex = workingIndex;
+ }
+
+ return true;
+ }
+ return false;
+ }
+
+ /**
+ * Move to the upper level
+ * It may cause to change the suboctree we are working on
+ * but we are using the same morton index >> 3
+ * @return true if succeed
+ */
+ bool moveUp() {
+ // It is on the top level?
+ if( this->currentLocalLevel ){
+ // No so simply go up
+ --this->currentLocalLevel;
+ this->currentLocalIndex >>= 3;
+ }
+ // Yes need to change suboctree
+ else if( this->current.tree->hasParent() ){
+ this->currentLocalIndex = this->current.tree->getIndexInParent();
+ this->current.tree = this->current.tree->getParent();
+ this->currentLocalLevel = this->current.tree->getSubOctreeHeight() - 1;
+ }
+ else{
+ return false;
+ }
+ return true;
+ }
+
+ /**
+ * Move down
+ * It may cause to change the suboctree we are working on
+ * We point on the first child found from left to right in the above
+ * level
+ * @return true if succeed
+ */
+ bool moveDown(){
+ if( !isAtLeafLevel() ){
+ // We are on the leaf of the current suboctree?
+ if(this->currentLocalLevel + 1 == this->current.tree->getSubOctreeHeight()){
+ // Yes change sub octree
+ this->current.tree = this->current.middleTree->leafs(this->currentLocalIndex);
+ this->currentLocalIndex = 0;
+ this->currentLocalLevel = 0;
+ }
+ // No simply go down
+ else{
+ ++this->currentLocalLevel;
+ this->currentLocalIndex <<= 3;
+ }
+ // Find the first allocated cell from left to right
+ while(!this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]){
+ ++this->currentLocalIndex;
+ }
+ return true;
+ }
+ return false;
+ }
+
+ /**
+ * To know if we are not on the root level
+ * @return true if we can move up
+ */
+ bool canProgressToUp() const {
+ return this->currentLocalLevel || this->current.tree->hasParent();
+ }
+
+ /**
+ * To know if we are not on the leafs level
+ * @return true if we can move down
+ */
+ bool canProgressToDown() const {
+ return !isAtLeafLevel();
+ }
+
+ /**
+ * To know if we are on the leafs level
+ * @return true if we are at the bottom of the tree
+ */
+ bool isAtLeafLevel() const {
+ return this->current.tree->isLeafPart() && this->currentLocalLevel + 1 == this->current.tree->getSubOctreeHeight();
+ }
+
+ /**
+ * To know the current level (not local but global)
+ * @return the level in the entire octree
+ */
+ int level() const {
+ return this->currentLocalLevel + this->current.tree->getSubOctreePosition();
+ }
+
+ /** Get the current pointed leaf
+ * @return current leaf element
+ */
+ LeafClass* getCurrentLeaf() const {
+ return this->current.leafTree->getLeaf(this->currentLocalIndex);
+ }
+
+ /** To access the current particles list
+ * You have to be at the leaf level to call this function!
+ * @return current element list
+ */
+ ContainerClass* getCurrentListSrc() const {
+ return this->current.leafTree->getLeafSrc(this->currentLocalIndex);
+ }
+
+ /** To access the current particles list
+ * You have to be at the leaf level to call this function!
+ * @return current element list
+ */
+ ContainerClass* getCurrentListTargets() const {
+ return this->current.leafTree->getLeafTargets(this->currentLocalIndex);
+ }
+
+ /** Get the current pointed cell
+ * @return current cell element
+ */
+ CellClass* getCurrentCell() const {
+ return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex];
+ }
+
+ /** Get the child of the current cell
+ * This function return an array of CellClass (array size = 8)
+ * User has to test each case to know if there is a cell
+ * @return the child array
+ */
+ CellClass** getCurrentChild() const {
+ // are we at the bottom of the suboctree
+ if(this->current.tree->getSubOctreeHeight() - 1 == this->currentLocalLevel ){
+ // then return first level of the suboctree under
+ return &this->current.middleTree->leafs(this->currentLocalIndex)->cellsAt(0)[0];
+ }
+ else{
+ // else simply return the array at the right position
+ return &this->current.tree->cellsAt(this->currentLocalLevel + 1)[this->currentLocalIndex << 3];
+ }
+ }
+
+ /** Get the part of array that contains all the pointers
+ *
+ */
+ CellClass** getCurrentBox() const {
+ return &this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex & ~7];
+ }
+
+ /** Get the Morton index of the current cell pointed by the iterator
+ * @return The global morton index
+ * <code>iter.getCurrentGlobalIndex();<br>
+ * // is equivalent to :<br>
+ * iter.getCurrentCell()->getMortonIndex();</code>
+ */
+ MortonIndex getCurrentGlobalIndex() const{
+ return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]->getMortonIndex();
+ }
+
+ /** To get the tree coordinate of the current working cell
+ *
+ */
+ const FTreeCoordinate& getCurrentGlobalCoordinate() const{
+ return this->current.tree->cellsAt(this->currentLocalLevel)[this->currentLocalIndex]->getCoordinate();
+ }
+
+ };
+
+ // To be able to access octree root & data
+ 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 find a missing cell or the right cell
+ */
+ CellClass* getCell(const MortonIndex inIndex, const int inLevel) const{
+ SubOctreeTypesConst workingTree;
+ workingTree.tree = &this->root;
+ const MortonIndex treeSubLeafMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
+
+ // Find the suboctree a the correct level
+ while(inLevel >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
+ // compute the leaf index
+ const MortonIndex fullIndex = inIndex >> (3 * (inLevel + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) ) );
+ // point to next suboctree
+ workingTree.tree = workingTree.middleTree->leafs(treeSubLeafMask & fullIndex);
+ if(!workingTree.tree) return 0;
+ }
+
+ // compute correct index in the array
+ const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (inLevel + 1 - workingTree.tree->getSubOctreePosition()) ));
+ return workingTree.tree->cellsAt(inLevel - workingTree.tree->getSubOctreePosition())[treeLeafMask & inIndex];
+ }
+
+
+ /** This function fill a array with the neighbors of a cell
+ * it does not put the brothers in the array (brothers are cells
+ * at the same level with the same parent) because they are of course
+ * direct neighbors.
+ * There is a maximum of 26 (3*3*3-1) direct neighbors
+ * // Take the neighbors != brothers
+ * CellClass* directNeighbors[26];
+ * const int nbDirectNeighbors = getNeighborsNoBrothers(directNeighbors,inIndex,inLevel);
+ * @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 getNeighborsNoBrothers(CellClass* inNeighbors[26], const MortonIndex inIndex, const int inLevel) const {
+ FTreeCoordinate center;
+ center.setPositionFromMorton(inIndex, inLevel);
+
+ const int boxLimite = FMath::pow2(inLevel);
+
+ int idxNeighbors = 0;
+
+ // We test all cells around
+ for(int idxX = -1 ; idxX <= 1 ; ++idxX){
+ if(!FMath::Between(center.getX() + idxX,0,boxLimite)) continue;
+
+ for(int idxY = -1 ; idxY <= 1 ; ++idxY){
+ if(!FMath::Between(center.getY() + idxY,0,boxLimite)) continue;
+
+ for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
+ if(!FMath::Between(center.getZ() + idxZ,0,boxLimite)) continue;
+
+ // 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(inLevel);
+ // if not a brother
+ if( mortonOther>>3 != inIndex>>3 ){
+ // get cell
+ CellClass* const cell = getCell(mortonOther, inLevel);
+ // add to list if not null
+ if(cell) inNeighbors[idxNeighbors++] = cell;
+ }
+ }
+ }
+ }
+ }
+
+ return idxNeighbors;
+ }
+
+
+ /** This function return an adresse of cell array from a morton index and a level
+ *
+ * @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)
+ *
+ */
+ CellClass** getCellPt(const MortonIndex inIndex, const int inLevel) const{
+ SubOctreeTypesConst workingTree;
+ workingTree.tree = this->root;
+
+ const MortonIndex treeMiddleMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
+
+ // Find the suboctree a the correct level
+ while(inLevel >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
+ // compute the leaf index
+ const MortonIndex fullIndex = inIndex >> 3 * (inLevel + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) );
+ // point to next suboctree
+ workingTree.tree = workingTree.middleTree->leafs(int(treeMiddleMask & fullIndex));
+ if(!workingTree.tree) return 0;
+ }
+
+ // Be sure there is a parent allocated
+ const int levelInTree = inLevel - workingTree.tree->getSubOctreePosition();
+ if( levelInTree && !workingTree.tree->cellsAt(levelInTree - 1)[~(~0x00LL << (3 * levelInTree )) & (inIndex>>3)]){
+ return 0;
+ }
+
+ // compute correct index in the array and return the @ in array
+ const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (levelInTree + 1 ) ));
+ return &workingTree.tree->cellsAt(levelInTree)[treeLeafMask & inIndex];
+ }
+
+
+ /** This function fill an array with the distant neighbors of a cell
+ * @param inNeighbors the array to store the elements
+ * @param inNeighborsIndex the array to store morton index of the neighbors
+ * @param inIndex the index of the element we want the neighbors
+ * @param inLevel the level of the element
+ * @return the number of neighbors
+ */
+ int getInteractionNeighbors(const CellClass* inNeighbors[343],
+ const FTreeCoordinate& workingCell,
+ const int inLevel) const{
+ // reset
+ memset(inNeighbors, 0, sizeof(CellClass*) * 343);
+
+ // Then take each child of the parent's neighbors if not in directNeighbors
+ // Father coordinate
+ const FTreeCoordinate parentCell(workingCell.getX()>>1,workingCell.getY()>>1,workingCell.getZ()>>1);
+
+ // Limite at parent level number of box (split by 2 by level)
+ const int boxLimite = FMath::pow2(inLevel-1);
+
+ int idxNeighbors = 0;
+ // We test all cells around
+ for(int idxX = -1 ; idxX <= 1 ; ++idxX){
+ if(!FMath::Between(parentCell.getX() + idxX,0,boxLimite)) continue;
+
+ for(int idxY = -1 ; idxY <= 1 ; ++idxY){
+ if(!FMath::Between(parentCell.getY() + idxY,0,boxLimite)) continue;
+
+ for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
+ if(!FMath::Between(parentCell.getZ() + idxZ,0,boxLimite)) continue;
+
+ // if we are not on the current cell
+ if( idxX || idxY || idxZ ){
+ const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
+ const MortonIndex mortonOtherParent = otherParent.getMortonIndex(inLevel-1) << 3;
+ // Get child
+ CellClass** const cells = getCellPt(mortonOtherParent, inLevel);
+
+ // If there is one or more child
+ if(cells){
+ // For each child
+ for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
+ if(cells[idxCousin]){
+ const int xdiff = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - workingCell.getX();
+ const int ydiff = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - workingCell.getY();
+ const int zdiff = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
+
+ // Test if it is a direct neighbor
+ if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+ // add to neighbors
+ inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
+ ++idxNeighbors;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return idxNeighbors;
+ }
+
+
+ /** This function fill an array with the distant neighbors of a cell
+ * it respects the periodic condition and will give the relative distance
+ * between the working cell and the neighbors
+ * @param inNeighbors the array to store the elements
+ * @param inRelativePosition the array to store the relative position of the neighbors
+ * @param workingCell the index of the element we want the neighbors
+ * @param inLevel the level of the element
+ * @return the number of neighbors
+ */
+ int getPeriodicInteractionNeighbors(const CellClass* inNeighbors[343],
+ const FTreeCoordinate& workingCell,
+ const int inLevel, const int inDirection) const{
+// std::cout << " Begin in getPeriodicInteractionNeighbors"<<std::endl;
+
+ // Then take each child of the parent's neighbors if not in directNeighbors
+ // Father coordinate
+ const FTreeCoordinate parentCell(workingCell.getX()>>1,workingCell.getY()>>1,workingCell.getZ()>>1);
+
+ // Limite at parent level number of box (split by 2 by level)
+ const int boxLimite = FMath::pow2(inLevel-1);
+
+ // This is not on a border we can use normal interaction list method
+ if( !(parentCell.getX() == 0 || parentCell.getY() == 0 || parentCell.getZ() == 0 ||
+ parentCell.getX() == boxLimite - 1 || parentCell.getY() == boxLimite - 1 || parentCell.getZ() == boxLimite - 1 ) ) {
+ return getInteractionNeighbors( inNeighbors, workingCell, inLevel);
+ }
+ else{
+ // reset
+
+ memset(inNeighbors, 0, sizeof(CellClass*) * 343);
+
+ const int startX = (TestPeriodicCondition(inDirection, DirMinusX) || parentCell.getX() != 0 ?-1:0);
+ const int endX = (TestPeriodicCondition(inDirection, DirPlusX) || parentCell.getX() != boxLimite - 1 ?1:0);
+ const int startY = (TestPeriodicCondition(inDirection, DirMinusY) || parentCell.getY() != 0 ?-1:0);
+ const int endY = (TestPeriodicCondition(inDirection, DirPlusY) || parentCell.getY() != boxLimite - 1 ?1:0);
+ const int startZ = (TestPeriodicCondition(inDirection, DirMinusZ) || parentCell.getZ() != 0 ?-1:0);
+ const int endZ = (TestPeriodicCondition(inDirection, DirPlusZ) || parentCell.getZ() != boxLimite - 1 ?1:0);
+// std::cout << " -- startX " << startX << " endX "<< endX<< std::endl ;
+// std::cout << " -- startY " << startY << " endX "<< endY<< std::endl ;
+// std::cout << " -- startZ " << startZ << " endX "<< endZ<< std::endl ;
+// std::cout << " boxLimite "<< boxLimite<<std::endl;
+ int idxNeighbors = 0;
+ // We test all cells around
+ for(int idxX = startX ; idxX <= endX ; ++idxX){
+ for(int idxY = startY ; idxY <= endY ; ++idxY){
+ for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
+ // if we are not on the current cell
+ if( idxX || idxY || idxZ ){
+
+ const FTreeCoordinate otherParent(parentCell.getX() + idxX,parentCell.getY() + idxY,parentCell.getZ() + idxZ);
+ FTreeCoordinate otherParentInBox(otherParent);
+ // periodic
+ if( otherParentInBox.getX() < 0 ){
+ otherParentInBox.setX( otherParentInBox.getX() + boxLimite );
+ }
+ else if( boxLimite <= otherParentInBox.getX() ){
+ otherParentInBox.setX( otherParentInBox.getX() - boxLimite );
+ }
+
+ if( otherParentInBox.getY() < 0 ){
+ otherParentInBox.setY( otherParentInBox.getY() + boxLimite );
+ }
+ else if( boxLimite <= otherParentInBox.getY() ){
+ otherParentInBox.setY( otherParentInBox.getY() - boxLimite );
+ }
+
+ if( otherParentInBox.getZ() < 0 ){
+ otherParentInBox.setZ( otherParentInBox.getZ() + boxLimite );
+ }
+ else if( boxLimite <= otherParentInBox.getZ() ){
+ otherParentInBox.setZ( otherParentInBox.getZ() - boxLimite );
+ }
+
+
+ const MortonIndex mortonOtherParent = otherParentInBox.getMortonIndex(inLevel-1) << 3;
+ // Get child
+ CellClass** const cells = getCellPt(mortonOtherParent, inLevel);
+// std::cout << " idx,idy,idz "<< idxX << " "<< idxY << " "<< idxZ
+// <<" mortonOtherParent"<< mortonOtherParent<< std::endl;
+
+ // If there is one or more child
+ if(cells){
+ // For each child
+ for(int idxCousin = 0 ; idxCousin < 8 ; ++idxCousin){
+ if(cells[idxCousin]){
+ const int xdiff = ((otherParent.getX()<<1) | ( (idxCousin>>2) & 1)) - workingCell.getX();
+ const int ydiff = ((otherParent.getY()<<1) | ( (idxCousin>>1) & 1)) - workingCell.getY();
+ const int zdiff = ((otherParent.getZ()<<1) | (idxCousin&1)) - workingCell.getZ();
+
+ // Test if it is a direct neighbor
+ if(FMath::Abs(xdiff) > 1 || FMath::Abs(ydiff) > 1 || FMath::Abs(zdiff) > 1){
+ // add to neighbors
+// std::cout << " Voisin numero "<< idxNeighbors
+// << " indexinTab "<< (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3
+// << " idxXousin " << idxCousin<< std::endl;
+ inNeighbors[ (((xdiff+3) * 7) + (ydiff+3)) * 7 + zdiff + 3] = cells[idxCousin];
+ ++idxNeighbors;
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+// std::cout << " End in getPeriodicInteractionNeighbors"<<std::endl;
+
+ return idxNeighbors;
+ }
+ }
+
+
+ /** 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)
+ *
+ */
+ ContainerClass* getLeafSrc(const MortonIndex inIndex){
+ SubOctreeTypes workingTree;
+ workingTree.tree = this->root;
+ const MortonIndex treeSubLeafMask = ~(~0x00LL << (3 * workingTree.tree->getSubOctreeHeight() ));
+
+ // Find the suboctree a the correct level
+ while(leafIndex >= workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) {
+ // compute the leaf index
+ const MortonIndex fullIndex = inIndex >> (3 * (leafIndex + 1 - (workingTree.tree->getSubOctreeHeight() + workingTree.tree->getSubOctreePosition()) ) );
+ // point to next suboctree
+ workingTree.tree = workingTree.middleTree->leafs(int(treeSubLeafMask & fullIndex));
+ if(!workingTree.tree) return 0;
+ }
+
+ // compute correct index in the array
+ const MortonIndex treeLeafMask = ~(~0x00LL << (3 * (leafIndex + 1 - workingTree.tree->getSubOctreePosition()) ));
+ return workingTree.leafTree->getLeafSrc(int(treeLeafMask & inIndex));
+ }
+
+ /** 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(ContainerClass* inNeighbors[27], const FTreeCoordinate& center, const int inLevel){
+ memset( inNeighbors, 0 , 27 * sizeof(ContainerClass*));
+ const int boxLimite = FMath::pow2(inLevel);
+
+ int idxNeighbors = 0;
+
+ // We test all cells around
+ for(int idxX = -1 ; idxX <= 1 ; ++idxX){
+ if(!FMath::Between(center.getX() + idxX,0,boxLimite)) continue;
+
+ for(int idxY = -1 ; idxY <= 1 ; ++idxY){
+ if(!FMath::Between(center.getY() + idxY,0,boxLimite)) continue;
+
+ for(int idxZ = -1 ; idxZ <= 1 ; ++idxZ){
+ if(!FMath::Between(center.getZ() + idxZ,0,boxLimite)) continue;
+
+ // 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(inLevel);
+ // get cell
+ ContainerClass* const leaf = getLeafSrc(mortonOther);
+ // add to list if not null
+ if(leaf){
+ inNeighbors[(((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1] = leaf;
+ ++idxNeighbors;
+ }
+ }
+ }
+ }
+ }
+
+ return idxNeighbors;
+ }
+
+
+ /** 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 getPeriodicLeafsNeighbors(ContainerClass* inNeighbors[27], FTreeCoordinate outOffsets[27], bool*const isPeriodic,
+ const FTreeCoordinate& center, const int inLevel, const int inDirection){
+
+ const int boxLimite = FMath::pow2(inLevel);
+// std::cout << " Begin in getPeriodicLeafsNeighbors"<<std::endl;
+
+ if( center.getX() != 0 && center.getY() != 0 && center.getZ() != 0 &&
+ center.getX() != boxLimite - 1 && center.getY() != boxLimite - 1 && center.getZ() != boxLimite - 1 ){
+ (*isPeriodic) = false;
+ return getLeafsNeighbors(inNeighbors, center, inLevel);
+ }
+
+ (*isPeriodic) = true;
+ memset(inNeighbors , 0 , 27 * sizeof(ContainerClass*));
+ int idxNeighbors = 0;
+
+ const int startX = (TestPeriodicCondition(inDirection, DirMinusX) || center.getX() != 0 ?-1:0);
+ const int endX = (TestPeriodicCondition(inDirection, DirPlusX) || center.getX() != boxLimite - 1 ?1:0);
+ const int startY = (TestPeriodicCondition(inDirection, DirMinusY) || center.getY() != 0 ?-1:0);
+ const int endY = (TestPeriodicCondition(inDirection, DirPlusY) || center.getY() != boxLimite - 1 ?1:0);
+ const int startZ = (TestPeriodicCondition(inDirection, DirMinusZ) || center.getZ() != 0 ?-1:0);
+ const int endZ = (TestPeriodicCondition(inDirection, DirPlusZ) || center.getZ() != boxLimite - 1 ?1:0);
+// std::cout << " -- startX " << startX << " endX "<< endX<< std::endl ;
+// std::cout << " -- startY " << startY << " endX "<< endY<< std::endl ;
+// std::cout << " -- startZ " << startZ << " endX "<< endZ<< std::endl ;
+// std::cout << " boxLimite "<< boxLimite<<std::endl;
+ int otherX,otherY,otherZ;
+ FTreeCoordinate other;
+
+ int xoffset = 0, yoffset = 0, zoffset = 0;
+ // We test all cells around
+ for(int idxX = startX ; idxX <= endX ; ++idxX){
+ otherX = center.getX() + idxX ; xoffset = 0 ;
+ if( otherX < 0 ){
+ otherX += boxLimite ;
+ xoffset = -1;
+ }
+ else if( boxLimite <= otherX ){
+ otherX -= boxLimite ;
+ xoffset = 1;
+ }
+ other.setX(otherX);
+ for(int idxY = startY ; idxY <= endY ; ++idxY){
+ otherY = center.getY() + idxY ;
+ yoffset = 0 ;
+ if( otherY < 0 ){
+ otherY += boxLimite ;
+ yoffset = -1;
+ }
+ else if( boxLimite <= otherY ){
+ otherY -= boxLimite ;
+ yoffset = 1;
+ }
+ other.setY(otherY);
+ for(int idxZ = startZ ; idxZ <= endZ ; ++idxZ){
+ zoffset = 0 ;
+ // if we are not on the current cell
+ if( idxX || idxY || idxZ ){ // !( idxX !=0 && idxY != 0 &&idxZ != 0 )
+// FTreeCoordinate other(center.getX() + idxX,center.getY() + idxY,center.getZ() + idxZ);
+ otherZ = center.getZ() + idxZ ;
+
+ if( otherZ < 0 ){
+ otherZ += boxLimite ;
+ zoffset = -1;
+ }
+ else if( boxLimite <= otherZ ){
+ otherZ -= boxLimite ;
+ zoffset = 1;
+ }
+ other.setZ(otherZ);
+
+// if( other.getX() < 0 ){
+// other.setX( other.getX() + boxLimite );
+// xoffset = -1;
+// }
+// else if( boxLimite <= other.getX() ){
+// other.setX( other.getX() - boxLimite );
+// xoffset = 1;
+// }
+// if( other.getY() < 0 ){
+// other.setY( other.getY() + boxLimite );
+// yoffset = -1;
+// }
+// else if( boxLimite <= other.getY() ){
+// other.setY( other.getY() - boxLimite );
+// yoffset = 1;
+// }
+// if( other.getZ() < 0 ){
+// other.setZ( other.getZ() + boxLimite );
+// zoffset = -1;
+// }
+// else if( boxLimite <= other.getZ() ){
+// other.setZ( other.getZ() - boxLimite );
+// zoffset = 1;
+// }
+
+
+ const MortonIndex mortonOther = other.getMortonIndex(inLevel);
+ //
+// std::cout << " idx,idy,idz "<< idxX << " "<< idxY << " "<< idxZ
+// <<" mortonOther "<< mortonOther<< " other "<< other<< std::endl;
+ // get cell
+ ContainerClass* const leaf = getLeafSrc(mortonOther);
+ // add to list if not null
+ if(leaf){
+ const int index = (((idxX + 1) * 3) + (idxY +1)) * 3 + idxZ + 1;
+ inNeighbors[index] = leaf;
+ outOffsets[index].setPosition(xoffset,yoffset,zoffset);
+// std::cout << " xoffset,yoffset,zoffset "<< xoffset << " "<< yoffset << " "<< zoffset <<" mortonOther " << mortonOther << std::endl<< std::endl;
+
+ ++idxNeighbors;
+ } //�if(leaf)
+ } // if( idxX || idxY || idxZ )
+ }
+ }
+ }
+// std::cout << " End in getPeriodicLeafsNeighbors " <<std::endl;
+
+ return idxNeighbors;
+ }
+
+
+
+ /////////////////////////////////////////////////////////
+ // Lambda function to apply to all member
+ /////////////////////////////////////////////////////////
+
+ /**
+ * @brief forEachLeaf iterate on the leaf and apply the function
+ * @param function
+ */
+ void forEachLeaf(std::function<void(LeafClass*)> function){
+ Iterator octreeIterator(this);
+ octreeIterator.gotoBottomLeft();
+
+ do{
+ function(octreeIterator.getCurrentLeaf());
+ } while(octreeIterator.moveRight());
+ }
+
+ /**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+ void forEachCell(std::function<void(CellClass*)> function){
+ Iterator octreeIterator(this);
+ octreeIterator.gotoBottomLeft();
+
+ Iterator avoidGoLeft(octreeIterator);
+
+ for(int idx = this->height-1 ; idx >= 1 ; --idx ){
+ do{
+ function(octreeIterator.getCurrentCell());
+ } while(octreeIterator.moveRight());
+ avoidGoLeft.moveUp();
+ octreeIterator = avoidGoLeft;
+ }
+ }
+
+ /**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+ void forEachCellWithLevel(std::function<void(CellClass*,const int)> function){
+ Iterator octreeIterator(this);
+ octreeIterator.gotoBottomLeft();
+
+ Iterator avoidGoLeft(octreeIterator);
+
+ for(int idx = this->height-1 ; idx >= 1 ; --idx ){
+ do{
+ function(octreeIterator.getCurrentCell(),idx);
+ } while(octreeIterator.moveRight());
+ avoidGoLeft.moveUp();
+ octreeIterator = avoidGoLeft;
+ }
+ }
+
+ /**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+ void forEachCellLeaf(std::function<void(CellClass*,LeafClass*)> function){
+ Iterator octreeIterator(this);
+ octreeIterator.gotoBottomLeft();
+
+ do{
+ function(octreeIterator.getCurrentCell(),octreeIterator.getCurrentLeaf());
+ } while(octreeIterator.moveRight());
+ }
};
#endif //FOCTREE_HPP
diff --git a/Src/Files/FEwalLoader.hpp b/Src/Files/FEwalLoader.hpp
index 989b00b69704be5412ed49c9aafb911f38d66980..0415d459ceebcb2f6940d5ac41620e9842f9d4ac 100755
--- a/Src/Files/FEwalLoader.hpp
+++ b/Src/Files/FEwalLoader.hpp
@@ -44,6 +44,8 @@ public:
enum Type{
OW,
HW,
+ Na,
+ Cl,
Undefined,
};
@@ -141,43 +143,56 @@ public:
-4406.48579000 6815.52906417 10340.2577024
*/
void fillParticle(FPoint* inPosition, FReal inForces[3], FReal* inPhysicalValue, int* inIndex){
- FReal x, y, z, fx, fy, fz, vx, vy, vz;
- int index;
- char type[2];
- std::string line;
- file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
-
- file.read(type, 2);
- file >> index;
- std::getline(file, line); // needed to skip the end of the line in non periodic case
- if ( levcfg == 0) {
- file >> x >> y >> z;
- }else if ( levcfg == 1) {
- file >> x >> y >> z;
- file >> vx >> vy >> vz;
- }else {
- file >> x >> y >> z;
- file >> vx >> vy >> vz;
- file >> fx >> fy >> fz;
- }
-
- inIndex[0] = index;
-
- inPosition->setPosition( x, y ,z);
- inForces[0] = fx;
- inForces[1] = fy;
- inForces[2] = fz;
-
- if( strncmp(type, "OW", 2) == 0){
- *inPhysicalValue = FReal(-0.82);
- *inIndex = OW;
- }
- else{
- *inPhysicalValue = FReal(-0.41);
- *inIndex = HW;
- }
+ FReal x, y, z, fx, fy, fz, vx, vy, vz;
+ int index;
+ char type[2];
+ std::string line;
+ file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+
+ file.read(type, 2);
+ file >> index;
+ std::getline(file, line); // needed to skip the end of the line in non periodic case
+ if ( levcfg == 0) {
+ file >> x >> y >> z;
+ }else if ( levcfg == 1) {
+ file >> x >> y >> z;
+ file >> vx >> vy >> vz;
+ }else {
+ file >> x >> y >> z;
+ file >> vx >> vy >> vz;
+ file >> fx >> fy >> fz;
+ }
+
+ inIndex[0] = index;
+
+ inPosition->setPosition( x, y ,z);
+ inForces[0] = fx;
+ inForces[1] = fy;
+ inForces[2] = fz;
+
+ if( strncmp(type, "OW", 2) == 0){
+ *inPhysicalValue = FReal(-0.82);
+ *inIndex = OW;
+ }
+ else if( strncmp(type, "HW", 2) == 0){
+ *inPhysicalValue = FReal(-0.41);
+ *inIndex = HW;
+ }
+ else if( strncmp(type, "Na", 2) == 0){
+ *inPhysicalValue = FReal(1.0);
+ *inIndex = Na;
+ }
+ else if( strncmp(type, "Cl", 2) == 0){
+ *inPhysicalValue = FReal(-1.0);
+ *inIndex = Cl;
+ }
+ else{
+ *inPhysicalValue = FReal(-0.41);
+ *inIndex = HW;
+ std::cerr << "Atom type not defined"<< std::endl;
+ exit(-1);
+ }
}
-
};
diff --git a/Src/Kernels/Chebyshev/FChebSymKernel.hpp b/Src/Kernels/Chebyshev/FChebSymKernel.hpp
index 3fc6bc3c3bd8373c0b6c597a34762b5568801cb1..8bd2deefb62ec8bafcf959c6bf971f91b7c0bb3c 100755
--- a/Src/Kernels/Chebyshev/FChebSymKernel.hpp
+++ b/Src/Kernels/Chebyshev/FChebSymKernel.hpp
@@ -123,6 +123,16 @@ public:
#endif
}
+ /**
+ * The constructor initializes all constant attributes and it reads the
+ * precomputed and compressed M2L operators from a binary file (an
+ * runtime_error is thrown if the required file is not valid).
+ */
+ FChebSymKernel(const int inTreeHeight,
+ const FReal inBoxWidth,
+ const FPoint& inBoxCenter) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter, FReal(1.0/FMath::pow(10,ORDER)))
+ {}
+
/** Copy constructor */