Remove some warnings, and split P2PR which is P2P only for one over R and P2P...

Remove some warnings, and split P2PR which is P2P only for one over R and P2P which can accept a matrix kernel (which can be one over R or anything else), both approaches support SSE and AVX but have not been tested.

Src/Files/FDlpolyLoader.hpp

@@ -300,7 +300,7 @@ public:
     */
     FDlpolyBinLoader(const char* const filename): file(fopen(filename, "rb")) {
         // test if open
-        if(this->file != NULL){
+        if(this->file != nullptr){
             energy = readValue<double>();
             double boxDim[3];
             boxWidth = readArray<double>(boxDim,3)[0];
@@ -325,7 +325,7 @@ public:
       * @return true if loader can work
       */
     bool isOpen() const{
-        return this->file != NULL;
+        return this->file != nullptr;
     }
 
 //    /**

Src/Kernels/P2P/FP2PClassic.hpp

 #ifndef FP2PCLASSIC_HPP
 #define FP2PCLASSIC_HPP
 
-#include "../../Utils/FGlobal.hpp"
-#include "../../Utils/FMath.hpp"
+namespace FP2P {
 
+/*
+   * FullMutual (generic version)
+   */
+template <class ContainerClass, typename MatrixKernelClass>
+inline void FullMutual(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
+                       const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
 
-namespace FP2P{
-
-  template <class ContainerClass>
-  inline void FullMutual(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
-			 const int limiteNeighbors){
 
     const int nbParticlesTargets = inTargets->getNbParticles();
     const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues();
@@ -22,79 +22,75 @@ namespace FP2P{
     FReal*const targetsPotentials = inTargets->getPotentials();
 
     for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
-      if( inNeighbors[idxNeighbors] ){
-	const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
-	const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
-	const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
-	const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
-	const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
-	FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX();
-	FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY();
-	FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ();
-	FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials();
-
-	for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
-	  for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
-	    FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
-	    FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
-	    FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
-
-	    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
-	    const FReal inv_distance = FMath::Sqrt(inv_square_distance);
-
-	    inv_square_distance *= inv_distance;
-	    inv_square_distance *= targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource];
-
-	    dx *= inv_square_distance;
-	    dy *= inv_square_distance;
-	    dz *= inv_square_distance;
-
-	    targetsForcesX[idxTarget] += dx;
-	    targetsForcesY[idxTarget] += dy;
-	    targetsForcesZ[idxTarget] += dz;
-	    targetsPotentials[idxTarget] += inv_distance * sourcesPhysicalValues[idxSource];
-
-	    sourcesForcesX[idxSource] -= dx;
-	    sourcesForcesY[idxSource] -= dy;
-	    sourcesForcesZ[idxSource] -= dz;
-	    sourcesPotentials[idxSource] += inv_distance * targetsPhysicalValues[idxTarget];
-	  }
-	}
-      }
+        for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
+            if( inNeighbors[idxNeighbors] ){
+                const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
+                const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
+                const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
+                const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
+                const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
+                FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX();
+                FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY();
+                FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ();
+                FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials();
+
+                for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
+
+                    // Compute kernel of interaction and its derivative
+                    const FPoint sourcePoint(sourcesX[idxSource],sourcesY[idxSource],sourcesZ[idxSource]);
+                    const FPoint targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
+                    FReal Kxy[1];
+                    FReal dKxy[3];
+                    MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
+                                                             targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
+                    FReal coef = (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
+
+                    targetsForcesX[idxTarget] += dKxy[0] * coef;
+                    targetsForcesY[idxTarget] += dKxy[1] * coef;
+                    targetsForcesZ[idxTarget] += dKxy[2] * coef;
+                    targetsPotentials[idxTarget] += ( Kxy[0] * sourcesPhysicalValues[idxSource] );
+
+                    sourcesForcesX[idxSource] -= dKxy[0] * coef;
+                    sourcesForcesY[idxSource] -= dKxy[1] * coef;
+                    sourcesForcesZ[idxSource] -= dKxy[2] * coef;
+                    sourcesPotentials[idxSource] += ( Kxy[0] * targetsPhysicalValues[idxTarget] );
+                }
+            }
+        }
     }
 
     for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
-      for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
-	FReal dx = targetsX[idxSource] - targetsX[idxTarget];
-	FReal dy = targetsY[idxSource] - targetsY[idxTarget];
-	FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
-
-	FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
-	const FReal inv_distance = FMath::Sqrt(inv_square_distance);
-
-	inv_square_distance *= inv_distance;
-	inv_square_distance *= targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource];
-
-	dx *= inv_square_distance;
-	dy *= inv_square_distance;
-	dz *= inv_square_distance;
-
-	targetsForcesX[idxTarget] += dx;
-	targetsForcesY[idxTarget] += dy;
-	targetsForcesZ[idxTarget] += dz;
-	targetsPotentials[idxTarget] += inv_distance * targetsPhysicalValues[idxSource];
-
-	targetsForcesX[idxSource] -= dx;
-	targetsForcesY[idxSource] -= dy;
-	targetsForcesZ[idxSource] -= dz;
-	targetsPotentials[idxSource] += inv_distance * targetsPhysicalValues[idxTarget];
-      }
+        for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
+
+            // Compute kernel of interaction...
+            const FPoint sourcePoint(targetsX[idxSource],targetsY[idxSource],targetsZ[idxSource]);
+            const FPoint targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
+            FReal Kxy[1];
+            FReal dKxy[3];
+            MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
+                                                     targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
+            FReal coef = (targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource]);
+
+            targetsForcesX[idxTarget] += dKxy[0] * coef;
+            targetsForcesY[idxTarget] += dKxy[1] * coef;
+            targetsForcesZ[idxTarget] += dKxy[2] * coef;
+            targetsPotentials[idxTarget] += ( Kxy[0] * targetsPhysicalValues[idxSource] );
+
+            targetsForcesX[idxSource] -= dKxy[0] * coef;
+            targetsForcesY[idxSource] -= dKxy[1] * coef;
+            targetsForcesZ[idxSource] -= dKxy[2] * coef;
+            targetsPotentials[idxSource] += ( Kxy[0] * targetsPhysicalValues[idxTarget] );
+        }
     }
-  }
+}
+
 
-template <class ContainerClass>
-  static void FullRemote(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
-			 const int limiteNeighbors){
+/**
+   * FullRemote (generic version)
+   */
+template <class ContainerClass, typename MatrixKernelClass>
+inline void FullRemote(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
+                       const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
     const int nbParticlesTargets = inTargets->getNbParticles();
     const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues();
     const FReal*const targetsX = inTargets->getPositions()[0];
@@ -106,39 +102,36 @@ template <class ContainerClass>
     FReal*const targetsPotentials = inTargets->getPotentials();
 
     for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
-      if( inNeighbors[idxNeighbors] ){
-	const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
-	const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
-	const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
-	const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
-	const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
-
-	for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
-
-	  for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
-	    FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
-	    FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
-	    FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
-
-	    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
-	    const FReal inv_distance = FMath::Sqrt(inv_square_distance);
-
-	    inv_square_distance *= inv_distance;
-	    inv_square_distance *= targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource];
-
-	    dx *= inv_square_distance;
-	    dy *= inv_square_distance;
-	    dz *= inv_square_distance;
-
-	    targetsForcesX[idxTarget] += dx;
-	    targetsForcesY[idxTarget] += dy;
-	    targetsForcesZ[idxTarget] += dz;
-	    targetsPotentials[idxTarget] += inv_distance * sourcesPhysicalValues[idxSource];
-	  }
-	}
-      }
+        for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
+            if( inNeighbors[idxNeighbors] ){
+                const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
+                const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
+                const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
+                const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
+                const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
+
+                for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
+
+                    // Compute kernel of interaction...
+                    const FPoint sourcePoint(sourcesX[idxSource],sourcesY[idxSource],sourcesZ[idxSource]);
+                    const FPoint targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
+                    FReal Kxy[1];
+                    FReal dKxy[3];
+                    MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
+                                                             targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
+                    FReal coef = (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
+
+                    targetsForcesX[idxTarget] += dKxy[0] * coef;
+                    targetsForcesY[idxTarget] += dKxy[1] * coef;
+                    targetsForcesZ[idxTarget] += dKxy[2] * coef;
+                    targetsPotentials[idxTarget] += ( Kxy[0] * sourcesPhysicalValues[idxSource] );
+                }
+            }
+        }
     }
-  }
+}
+
 
 }
-#endif //FP2PCLASSIC_HPP
+
+#endif // FP2PCLASSIC_HPP

Src/Kernels/P2P/FP2PR.hpp

+#ifndef FP2PR_HPP
+#define FP2PR_HPP
+
+#include "Utils/FGlobal.hpp"
+#include "Utils/FMath.hpp"
+
+
+/**
+ * @brief The FP2PR namespace
+ */
+namespace FP2PR{
+inline void MutualParticles(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue,
+                            FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
+                            const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
+                            FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential
+                            ){
+    FReal dx = sourceX - targetX;
+    FReal dy = sourceY - targetY;
+    FReal dz = sourceZ - targetZ;
+
+    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
+    FReal inv_distance = FMath::Sqrt(inv_square_distance);
+
+    inv_square_distance *= inv_distance;
+    inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
+
+    dx *= inv_square_distance;
+    dy *= inv_square_distance;
+    dz *= inv_square_distance;
+
+    *targetForceX += dx;
+    *targetForceY += dy;
+    *targetForceZ += dz;
+    *targetPotential += ( inv_distance * sourcePhysicalValue );
+
+    *sourceForceX -= dx;
+    *sourceForceY -= dy;
+    *sourceForceZ -= dz;
+    *sourcePotential += ( inv_distance * targetPhysicalValue );
+}
+
+inline void NonMutualParticles(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue,
+                               const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
+                               FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential){
+    FReal dx = sourceX - targetX;
+    FReal dy = sourceY - targetY;
+    FReal dz = sourceZ - targetZ;
+
+    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
+    FReal inv_distance = FMath::Sqrt(inv_square_distance);
+
+    inv_square_distance *= inv_distance;
+    inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
+
+    dx *= inv_square_distance;
+    dy *= inv_square_distance;
+    dz *= inv_square_distance;
+
+    *targetForceX += dx;
+    *targetForceY += dy;
+    *targetForceZ += dz;
+    *targetPotential += ( inv_distance * sourcePhysicalValue );
+}
+
+}
+
+#ifdef ScalFMM_USE_SSE
+#include "FP2PRSse.hpp"
+#elif defined(ScalFMM_USE_AVX)
+#include "FP2PRAvx.hpp"
+#else
+#include "FP2PRClassic.hpp"
+#endif //Includes
+
+#endif // FP2PR_HPP

Src/Kernels/P2P/FP2PAvx.hpp → Src/Kernels/P2P/FP2PRAvx.hpp

-#ifndef FP2PAVX_HPP
-#define FP2PAVX_HPP
+#ifndef FP2PRAVX_HPP
+#define FP2PRAVX_HPP
 
 #include "../../Utils/FGlobal.hpp"
 #include "../../Utils/FMath.hpp"
@@ -7,7 +7,7 @@
 #include "../../Utils/FAvx.hpp"
 
 
-namespace FP2P{
+namespace FP2PR{
 
 #ifdef ScalFMM_USE_DOUBLE_PRECISION
     template <class ContainerClass>
@@ -513,4 +513,4 @@ namespace FP2P{
 
 #endif
 }
-#endif //FP2PAVX_HPP
+#endif //FP2PRAVX_HPP

Src/Kernels/P2P/FP2PRClassic.hpp

+#ifndef FP2PRCLASSIC_HPP
+#define FP2PRCLASSIC_HPP
+
+#include "../../Utils/FGlobal.hpp"
+#include "../../Utils/FMath.hpp"
+
+
+namespace FP2PR{
+
+template <class ContainerClass>
+inline void FullMutual(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
+                       const int limiteNeighbors){
+
+    const int nbParticlesTargets = inTargets->getNbParticles();
+    const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues();
+    const FReal*const targetsX = inTargets->getPositions()[0];
+    const FReal*const targetsY = inTargets->getPositions()[1];
+    const FReal*const targetsZ = inTargets->getPositions()[2];
+    FReal*const targetsForcesX = inTargets->getForcesX();
+    FReal*const targetsForcesY = inTargets->getForcesY();
+    FReal*const targetsForcesZ = inTargets->getForcesZ();
+    FReal*const targetsPotentials = inTargets->getPotentials();
+
+    for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
+        if( inNeighbors[idxNeighbors] ){
+            const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
+            const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
+            const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
+            const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
+            const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
+            FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX();
+            FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY();
+            FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ();
+            FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials();
+
+            for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
+                for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
+                    FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
+                    FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
+                    FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
+
+                    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
+                    const FReal inv_distance = FMath::Sqrt(inv_square_distance);
+
+                    inv_square_distance *= inv_distance;
+                    inv_square_distance *= targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource];
+
+                    dx *= inv_square_distance;
+                    dy *= inv_square_distance;
+                    dz *= inv_square_distance;
+
+                    targetsForcesX[idxTarget] += dx;
+                    targetsForcesY[idxTarget] += dy;
+                    targetsForcesZ[idxTarget] += dz;
+                    targetsPotentials[idxTarget] += inv_distance * sourcesPhysicalValues[idxSource];
+
+                    sourcesForcesX[idxSource] -= dx;
+                    sourcesForcesY[idxSource] -= dy;
+                    sourcesForcesZ[idxSource] -= dz;
+                    sourcesPotentials[idxSource] += inv_distance * targetsPhysicalValues[idxTarget];
+                }
+            }
+        }
+    }
+
+    for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
+        for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
+            FReal dx = targetsX[idxSource] - targetsX[idxTarget];
+            FReal dy = targetsY[idxSource] - targetsY[idxTarget];
+            FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
+
+            FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
+            const FReal inv_distance = FMath::Sqrt(inv_square_distance);
+
+            inv_square_distance *= inv_distance;
+            inv_square_distance *= targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource];
+
+            dx *= inv_square_distance;
+            dy *= inv_square_distance;
+            dz *= inv_square_distance;
+
+            targetsForcesX[idxTarget] += dx;
+            targetsForcesY[idxTarget] += dy;
+            targetsForcesZ[idxTarget] += dz;
+            targetsPotentials[idxTarget] += inv_distance * targetsPhysicalValues[idxSource];
+
+            targetsForcesX[idxSource] -= dx;
+            targetsForcesY[idxSource] -= dy;
+            targetsForcesZ[idxSource] -= dz;
+            targetsPotentials[idxSource] += inv_distance * targetsPhysicalValues[idxTarget];
+        }
+    }
+}
+
+template <class ContainerClass>
+static void FullRemote(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
+                       const int limiteNeighbors){
+    const int nbParticlesTargets = inTargets->getNbParticles();
+    const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues();
+    const FReal*const targetsX = inTargets->getPositions()[0];
+    const FReal*const targetsY = inTargets->getPositions()[1];
+    const FReal*const targetsZ = inTargets->getPositions()[2];
+    FReal*const targetsForcesX = inTargets->getForcesX();
+    FReal*const targetsForcesY = inTargets->getForcesY();
+    FReal*const targetsForcesZ = inTargets->getForcesZ();
+    FReal*const targetsPotentials = inTargets->getPotentials();
+
+    for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
+        if( inNeighbors[idxNeighbors] ){
+            const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
+            const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
+            const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
+            const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
+            const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
+
+            for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
+
+                for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
+                    FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
+                    FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
+                    FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
+
+                    FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
+                    const FReal inv_distance = FMath::Sqrt(inv_square_distance);
+
+                    inv_square_distance *= inv_distance;
+                    inv_square_distance *= targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource];
+
+                    dx *= inv_square_distance;
+                    dy *= inv_square_distance;
+                    dz *= inv_square_distance;
+
+                    targetsForcesX[idxTarget] += dx;
+                    targetsForcesY[idxTarget] += dy;
+                    targetsForcesZ[idxTarget] += dz;
+                    targetsPotentials[idxTarget] += inv_distance * sourcesPhysicalValues[idxSource];
+                }
+            }
+        }
+    }
+}
+
+}
+#endif //FP2PRCLASSIC_HPP

Src/Kernels/Rotation/FRotationKernel.hpp

@@ -23,7 +23,7 @@
 #include "Utils/FSpherical.hpp"
 #include "Utils/FMath.hpp"
 
-#include "../P2P/FP2P.hpp"
+#include "../P2P/FP2PR.hpp"
 
 /** This is a recursion to get the minimal size of the matrix dlmk
   */
@@ -1310,7 +1310,7 @@ public:
     void P2P(const FTreeCoordinate& /*inPosition*/,
                      ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                      ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        FP2P::FullMutual(inTargets,inNeighbors,14);
+        FP2PR::FullMutual(inTargets,inNeighbors,14);
     }
 
 
@@ -1318,7 +1318,7 @@ public:
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        FP2P::FullRemote(inTargets,inNeighbors,27);
+        FP2PR::FullRemote(inTargets,inNeighbors,27);
     }
 };
 

Src/Kernels/Rotation/FRotationOriginalKernel.hpp

@@ -22,7 +22,7 @@
 #include "../../Utils/FMemUtils.hpp"
 #include "../../Utils/FSpherical.hpp"
 
-#include "../P2P/FP2P.hpp"
+#include "../P2P/FP2PR.hpp"
 
 /**
 * @author Berenger Bramas (berenger.bramas@inria.fr)
@@ -800,7 +800,7 @@ public:
     void P2P(const FTreeCoordinate& /*inPosition*/,
                      ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                      ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        FP2P::FullMutual(inTargets,inNeighbors,14);
+        FP2PR::FullMutual(inTargets,inNeighbors,14);
     }
 
 
@@ -808,7 +808,7 @@ public:
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        FP2P::FullRemote(inTargets,inNeighbors,27);
+        FP2PR::FullRemote(inTargets,inNeighbors,27);
     }
 };
 

Src/Kernels/Spherical/FAbstractSphericalKernel.hpp

@@ -28,7 +28,7 @@
 
 #include "../../Containers/FTreeCoordinate.hpp"
 
-#include "../P2P/FP2P.hpp"
+#include "../P2P/FP2PR.hpp"
 
 #include "FHarmonic.hpp"
 
@@ -222,7 +222,7 @@ public:
     void P2P(const FTreeCoordinate& inLeafPosition,
                   ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                   ContainerClass* const directNeighborsParticles[27], const int /*size*/){
-        FP2P::FullMutual(targets,directNeighborsParticles,14);
+        FP2PR::FullMutual(targets,directNeighborsParticles,14);
     }
 
     /** This P2P has to be used when target != sources
@@ -236,7 +236,7 @@ public:
     void P2PRemote(const FTreeCoordinate& ,
                   ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,
                   ContainerClass* const directNeighborsParticles[27], const int /*size*/){
-        FP2P::FullRemote(targets,directNeighborsParticles,27);
+        FP2PR::FullRemote(targets,directNeighborsParticles,27);
     }
 
 private:

Src/Kernels/Taylor/FTaylorKernel.hpp

@@ -20,7 +20,7 @@