Implemented dense tensorial Chebyshev kernel + add test, fixed some failing...

Implemented dense tensorial Chebyshev kernel + add test, fixed some failing utests by comparing RELATIVE errors instead of absolute.

Src/Kernels/Chebyshev/FAbstractChebKernel.hpp

@@ -45,7 +45,7 @@ class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass>
 {
 protected:
   enum {nnodes = TensorTraits<ORDER>::nnodes};
-  typedef FChebInterpolator<ORDER> InterpolatorClass;
+  typedef FChebInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
 
   /// Needed for P2M, M2M, L2L and L2P operators
   const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
@@ -59,6 +59,8 @@ protected:
   const FReal BoxWidth;    
   /// Width of a leaf cell box 
   const FReal BoxWidthLeaf;
+  /// Parameter to pass to matrix kernel (material specific or anything)
+  const double MatParam;
 
   /**
    * Compute center of leaf cell from its tree coordinate.
@@ -79,14 +81,16 @@ public:
    * runtime_error is thrown if the required file is not valid).
    */
   FAbstractChebKernel(const int inTreeHeight,
-		      const FReal inBoxWidth,
-		      const FPoint& inBoxCenter)
+                      const FReal inBoxWidth,
+                      const FPoint& inBoxCenter,
+                      const double inMatParam = 0.0)
     : Interpolator(new InterpolatorClass()),
-      MatrixKernel(new MatrixKernelClass()),
+      MatrixKernel(new MatrixKernelClass(inMatParam)),
       TreeHeight(inTreeHeight),
       BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
       BoxWidth(inBoxWidth),
-      BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
+      BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
+      MatParam(inMatParam)
   {
     /* empty */
   }

Src/Kernels/Chebyshev/FChebTensorialKernel.hpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.  
+// 
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info". 
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+#ifndef FCHEBTENSORIALKERNEL_HPP
+#define FCHEBTENSORIALKERNEL_HPP
+
+#include "../../Utils/FGlobal.hpp"
+#include "../../Utils/FTrace.hpp"
+#include "../../Utils/FSmartPointer.hpp"
+
+#include "./FAbstractChebKernel.hpp"
+//#include "./FChebM2LHandler.hpp"
+#include "./FChebTensorialM2LHandler.hpp" //PB: temporary version
+
+class FTreeCoordinate;
+
+/**
+ * @author Matthias Messner(matthias.messner@inria.fr)
+ * @class FChebTensorialKernel
+ * @brief
+ * Please read the license
+ *
+ * This kernels implement the Chebyshev interpolation based FMM operators. It
+ * implements all interfaces (P2P, P2M, M2M, M2L, L2L, L2P) which are required by
+ * the FFmmAlgorithm and FFmmAlgorithmThread.
+ *
+ * @tparam CellClass Type of cell
+ * @tparam ContainerClass Type of container to store particles
+ * @tparam MatrixKernelClass Type of matrix kernel function
+ * @tparam ORDER Chebyshev interpolation order
+ */
+template < class CellClass,	class ContainerClass,	class MatrixKernelClass, int ORDER, int NVALS = 1>
+class FChebTensorialKernel
+    : public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
+{
+  enum {nRhs = MatrixKernelClass::NRHS,
+        nLhs = MatrixKernelClass::NLHS};
+
+protected://PB: for OptiDis
+
+	// private types
+	typedef FChebTensorialM2LHandler<ORDER,MatrixKernelClass,MatrixKernelClass::Type> M2LHandlerClass;
+
+	// using from 
+    typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
+	AbstractBaseClass;
+
+	/// Needed for M2L operator
+	FSmartPointer<  M2LHandlerClass,FSmartPointerMemory> M2LHandler;
+
+public:
+	/**
+	 * 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).
+	 */
+	FChebTensorialKernel(const int inTreeHeight,
+                       const FReal inBoxWidth,
+                       const FPoint& inBoxCenter,
+                       const FReal Epsilon,
+                       const double inMatParam = 0.0)
+        : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inMatParam),
+	  M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
+                                   inTreeHeight,
+                                   inBoxWidth,
+                                   Epsilon))
+	{ }
+
+
+	void P2M(CellClass* const LeafCell,
+					 const ContainerClass* const SourceParticles)
+	{
+    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
+
+    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+
+      // 1) apply Sy
+      AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
+                                                LeafCell->getMultipole(idxV*nRhs), SourceParticles);
+
+      for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
+        // update multipole index
+        int idxMul = idxV*nRhs + idxRhs;
+
+        // 2) apply B (PB: Tensorial version is just a basic copy)
+        M2LHandler->applyB(LeafCell->getMultipole(idxMul), LeafCell->getMultipole(idxMul) + AbstractBaseClass::nnodes);
+
+      }
+    }
+	}
+
+
+	void M2M(CellClass* const FRestrict ParentCell,
+					 const CellClass*const FRestrict *const FRestrict ChildCells,
+           const int /*TreeLevel*/)
+	{
+    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+      for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
+        // update multipole index
+        int idxMul = idxV*nRhs + idxRhs;
+
+        // 1) apply Sy
+        FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxMul));
+        for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
+          if (ChildCells[ChildIndex]){
+            AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxMul),
+                                                      ParentCell->getMultipole(idxMul));
+          }
+        }
+
+        // 2) apply B (PB: Tensorial version is just a basic copy)
+        M2LHandler->applyB(ParentCell->getMultipole(idxMul), ParentCell->getMultipole(idxMul) + AbstractBaseClass::nnodes);
+
+      }
+    }
+	}
+
+	void M2L(CellClass* const FRestrict TargetCell,
+					 const CellClass* SourceCells[343],
+           const int /*NumSourceCells*/,
+					 const int TreeLevel)
+	{
+    const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
+    const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth));
+
+    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+      for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
+        // update local index
+        int idxLoc = idxV*nLhs + idxLhs;
+
+        FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxLoc) + AbstractBaseClass::nnodes;
+      
+        for (int idxRhs=0; idxRhs < nRhs; ++idxRhs){
+          // update multipole index
+          int idxMul = idxV*nRhs + idxRhs;
+          // update kernel index such that: x_i = K_{ij}y_j 
+          int idxK = idxLhs*nRhs + idxRhs;
+          // get index in matrix kernel
+          unsigned int d 
+            = AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxK);
+
+          for (int idx=0; idx<343; ++idx){
+            if (SourceCells[idx]){
+              M2LHandler->applyC(idx, TreeLevel, scale, d,
+                                 SourceCells[idx]->getMultipole(idxMul) + AbstractBaseClass::nnodes,
+                                 CompressedLocalExpansion);
+            }
+          }
+        }// NRHS
+      }// NLHS
+    }// NVALS
+
+  }
+
+	void L2L(const CellClass* const FRestrict ParentCell,
+					 CellClass* FRestrict *const FRestrict ChildCells,
+           const int /*TreeLevel*/)
+	{
+    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+      for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
+        int idxLoc = idxV*nLhs + idxLhs;
+
+        // 1) apply U (PB: Tensorial version is just a basic copy)
+        M2LHandler->applyU(ParentCell->getLocal(idxLoc) + AbstractBaseClass::nnodes,
+                           const_cast<CellClass*>(ParentCell)->getLocal(idxLoc));
+        // 2) apply Sx
+        for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
+          if (ChildCells[ChildIndex]){
+            AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxLoc), ChildCells[ChildIndex]->getLocal(idxLoc));
+          }
+        }
+      }//NLHS
+    }// NVALS
+  }
+
+
+	void L2P(const CellClass* const LeafCell,
+					 ContainerClass* const TargetParticles)
+	{
+    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
+
+    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+      for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
+        int idxLoc = idxV*nLhs + idxLhs;
+
+        // 1) apply U (PB: Tensorial version is just a basic copy)
+        M2LHandler->applyU(LeafCell->getLocal(idxLoc) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxLoc));
+
+      }
+
+//      // 2.a) apply Sx
+//      AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
+//      																					AbstractBaseClass::BoxWidthLeaf,
+//      																					LeafCell->getLocal(idxV*nLhs),
+//      																					TargetParticles);
+//      // 2.b) apply Px (grad Sx)
+//      AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
+//      																									AbstractBaseClass::BoxWidthLeaf,
+//      																									LeafCell->getLocal(idxV*nLhs),
+//      																									TargetParticles);
+
+      // 2.c) apply Sx and Px (grad Sx)
+      AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
+                                                     LeafCell->getLocal(idxV*nLhs), TargetParticles);
+    }
+	}
+
+    void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
+                     ContainerClass* const FRestrict TargetParticles,
+                     const ContainerClass* const FRestrict /*SourceParticles*/,
+                     ContainerClass* const NeighborSourceParticles[27],
+                     const int /* size */)
+    {
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
+    }
+
+
+    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
+                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
+                   ContainerClass* const inNeighbors[27], const int /*inSize*/){
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
+    }
+
+};
+
+
+#endif //FCHEBTENSORIALKERNELS_HPP
+
+// [--END--]

Src/Kernels/Interpolation/FInterpMatrixKernel.hpp

@@ -89,6 +89,11 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
                                    xy.getZ()*xy.getZ());
   }
 
+  void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
+  {
+    block[0]=this->evaluate(x,y);
+  }
+
   FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
   {
     const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
@@ -210,8 +215,8 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R;
   static const unsigned int DIM = 6; //PB: dimension of kernel
   static const unsigned int NIDX = 2; //PB: number of indices 
-  static constexpr unsigned int indexTab[12]={0,0,0,1,1,2,
-                                   0,1,2,1,2,2};
+  /*static */const/*expr*/ unsigned int indexTab[12]={0,0,0,1,1,2,
+                                                      0,1,2,1,2,2};
   static const unsigned int NRHS = 3;
   static const unsigned int NLHS = 3;
 
@@ -280,8 +285,8 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJ;
   static const unsigned int DIM = 6; //PB: dimension of kernel
   static const unsigned int NIDX = 2; //PB: number of indices 
-  static constexpr unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,2,
-                                                    0,1,2,1,2,2};
+  /*static */const/*expr*/ unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,2,
+                                                            0,1,2,1,2,2};
   static const unsigned int NRHS = 3;
   static const unsigned int NLHS = 3;
 
@@ -364,7 +369,7 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJK;
   static const  unsigned int DIM = 10; //PB: dimension of kernel
   static const unsigned int NIDX = 3; //PB: number of indices 
-  static constexpr unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,1,2,2,2,0,
+  /*static */const/*expr*/ unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,1,2,2,2,0,
                                          0,1,2,0,1,2,0,1,2,1,
                                          0,1,2,0,1,2,0,1,2,2};
   static const unsigned int NRHS = 3;

Tests/Kernels/testChebTensorialAlgorithm.cpp

@@ -15,7 +15,7 @@
 // ===================================================================================
 
 // ==== CMAKE =====
-// @FUSE_BLAS
+// @FUSE_FFT
 // ================
 
 #include <iostream>
@@ -30,7 +30,7 @@
 
 #include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
 #include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
-#include "../../Src/Kernels/Chebyshev/FChebSymTensorialKernel.hpp"
+#include "../../Src/Kernels/Chebyshev/FChebTensorialKernel.hpp"
 
 #include "../../Src/Components/FSimpleLeaf.hpp"
 #include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
@@ -45,7 +45,7 @@
 #include "../../Src/Core/FFmmAlgorithmThread.hpp"
 
 /**
- * This program runs the FMM Algorithm with the Uniform kernel and compares the results with a direct computation.
+ * This program runs the FMM Algorithm with the Chebyshev kernel and compares the results with a direct computation.
  */
 
 // Simply create particles and try the kernels
@@ -66,14 +66,13 @@ int main(int argc, char* argv[])
     // init timer
     FTic time;
 
+
   // typedefs
-//  typedef FInterpMatrixKernelLJ MatrixKernelClass;
-//    typedef FInterpMatrixKernelR MatrixKernelClass;
-//    typedef FInterpMatrixKernel_R_IJ MatrixKernelClass; // not working with Non-Symmetric variant ! Because of UCB decomposition.
-                                                          // and not working either with Symmetric variant because not symmetric... 
-  typedef FInterpMatrixKernel_IOR MatrixKernelClass;
+  typedef FInterpMatrixKernel_R_IJ MatrixKernelClass;
+//  typedef FInterpMatrixKernel_IOR MatrixKernelClass;
+//  typedef FInterpMatrixKernelR MatrixKernelClass;
 
-//  const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
+  const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
   const unsigned int NRHS = MatrixKernelClass::NRHS;
   const unsigned int NLHS = MatrixKernelClass::NLHS;
 
@@ -97,7 +96,7 @@ int main(int argc, char* argv[])
     // get copy
     particles[idxPart].position       = position;
     for(unsigned idxRhs = 0; idxRhs<NRHS;++idxRhs)
-      particles[idxPart].physicalValue[idxRhs]  = physicalValue;
+      particles[idxPart].physicalValue[idxRhs]  = physicalValue; // copy same physical value in each component
     for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
       particles[idxPart].potential[idxLhs]      = 0.0;
       particles[idxPart].forces[0][idxLhs]      = 0.0;
@@ -114,22 +113,35 @@ int main(int argc, char* argv[])
     {
       for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
         for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
-//          FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
-//                                particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
-//                                &particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
-//                                &particles[idxTarget].forces[2], &particles[idxTarget].potential,
-//                                particles[idxOther].position.getX(), particles[idxOther].position.getY(),
-//                                particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
-//                                &particles[idxOther].forces[0], &particles[idxOther].forces[1],
-//                                &particles[idxOther].forces[2], &particles[idxOther].potential);
-          FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
-                                   particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
-                                   particles[idxTarget].forces[0], particles[idxTarget].forces[1],
-                                   particles[idxTarget].forces[2], particles[idxTarget].potential,
-                                   particles[idxOther].position.getX(), particles[idxOther].position.getY(),
-                                   particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
-                                   particles[idxOther].forces[0], particles[idxOther].forces[1],
-                                   particles[idxOther].forces[2], particles[idxOther].potential);
+          if(MK_ID == R_IJ)
+            FP2P::MutualParticlesRIJ(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
+                                     particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
+                                     particles[idxTarget].forces[0], particles[idxTarget].forces[1],
+                                     particles[idxTarget].forces[2], particles[idxTarget].potential,
+                                     particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+                                     particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
+                                     particles[idxOther].forces[0], particles[idxOther].forces[1],
+                                     particles[idxOther].forces[2], particles[idxOther].potential);
+          else if(MK_ID == ID_OVER_R)
+            FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
+                                     particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
+                                     particles[idxTarget].forces[0], particles[idxTarget].forces[1],
+                                     particles[idxTarget].forces[2], particles[idxTarget].potential,
+                                     particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+                                     particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
+                                     particles[idxOther].forces[0], particles[idxOther].forces[1],
+                                     particles[idxOther].forces[2], particles[idxOther].potential);
+          else if(MK_ID == ONE_OVER_R)
+            FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
+                                  particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue[0],
+                                  particles[idxTarget].forces[0], particles[idxTarget].forces[1],
+                                  particles[idxTarget].forces[2], particles[idxTarget].potential,
+                                  particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+                                  particles[idxOther].position.getZ(), particles[idxOther].physicalValue[0],
+                                  particles[idxOther].forces[0], particles[idxOther].forces[1],
+                                  particles[idxOther].forces[2], particles[idxOther].potential);
+          else 
+            std::runtime_error("Provide a valid matrix kernel!");
         }
       }
     }
@@ -141,19 +153,18 @@ int main(int argc, char* argv[])
 
   ////////////////////////////////////////////////////////////////////
 
-  {	// begin Chebyshev kernel
+  {	// begin Lagrange kernel
 
     // accuracy
-    const unsigned int ORDER = 4;
-    const FReal epsilon = FReal(1e-8);
+    const unsigned int ORDER = 6; //PB: Beware order=9 exceed memory (even 8 since compute _C then store in C) 
+    const FReal epsilon = FReal(1e-7);
 
-    // typedefs
     typedef FP2PParticleContainerIndexed<NRHS,NLHS> ContainerClass;
 
     typedef FSimpleLeaf< ContainerClass >  LeafClass;
     typedef FChebCell<ORDER,NRHS,NLHS> CellClass;
     typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
-    typedef FChebSymTensorialKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; // PB: tensorial kernel needs to be symmetric !!
+    typedef FChebTensorialKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
     typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
     //  typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
 
@@ -170,7 +181,13 @@ int main(int argc, char* argv[])
       for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
         // put in tree
         // PB: here we have to know the number of NRHS...
-        tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
+        if(NRHS==1)
+          tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue[0]);
+        else if(NRHS==3)
+          tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
+        else 
+          std::runtime_error("Insert correct number of physical values!");
+
       }
 
       time.tac();
@@ -202,6 +219,7 @@ int main(int argc, char* argv[])
 
         tree.forEachLeaf([&](LeafClass* leaf){
             for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
+
               const FReal*const physVals = leaf->getTargets()->getPhysicalValues(idxLhs);
               const FReal*const potentials = leaf->getTargets()->getPotentials(idxLhs);
               const FReal*const forcesX = leaf->getTargets()->getForcesX(idxLhs);
@@ -216,14 +234,14 @@ int main(int argc, char* argv[])
                 //PB: store potential in array[nbParticles]
                 checkPhysVal[indexPartOrig][idxLhs]=physVals[idxPart];              
                 checkPotential[indexPartOrig][idxLhs]=potentials[idxPart];              
-                checkfx[indexPartOrig][idxLhs]=forcesX[idxPart]; 
+                checkfx[indexPartOrig][idxLhs]=forcesX[idxPart];              
 
                 potentialDiff[idxLhs].add(particles[indexPartOrig].potential[idxLhs],potentials[idxPart]);
                 fx[idxLhs].add(particles[indexPartOrig].forces[0][idxLhs],forcesX[idxPart]);
                 fy[idxLhs].add(particles[indexPartOrig].forces[1][idxLhs],forcesY[idxPart]);
                 fz[idxLhs].add(particles[indexPartOrig].forces[2][idxLhs],forcesZ[idxPart]);
               }
-            }
+            }// NLHS
           });
       }
 
@@ -236,22 +254,20 @@ int main(int argc, char* argv[])
 //        }
 //      std::cout << std::endl;
 
-
       // Print for information
-      std::cout << "\nTest Chebyshev Tensorial returns wrong results! ChebInterpolator still needs to be modified like UnifInterpolator." << std::endl;
-//      std::cout << "\nAbsolute errors: " << std::endl;
-//      std::cout << "Potential: ";
-//      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << potentialDiff[idxLhs] << ", " ;
-//      std::cout << std::endl;
-//      std::cout << "Fx: "; 
-//      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fx[idxLhs] << ", " ;
-//      std::cout  << std::endl;
-//      std::cout << "Fy: "; 
-//      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fy[idxLhs] << ", " ;
-//      std::cout  << std::endl;
-//      std::cout << "Fz: "; 
-//      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fz[idxLhs] << ", " ;
-//      std::cout << std::endl;
+      std::cout << "\nAbsolute errors: " << std::endl;
+      std::cout << "Potential: ";
+      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << potentialDiff[idxLhs] << ", " ;
+      std::cout << std::endl;
+      std::cout << "Fx: "; 
+      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fx[idxLhs] << ", " ;
+      std::cout  << std::endl;
+      std::cout << "Fy: "; 
+      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fy[idxLhs] << ", " ;
+      std::cout  << std::endl;
+      std::cout << "Fz: "; 
+      for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fz[idxLhs] << ", " ;
+      std::cout << std::endl;
 
     } // -----------------------------------------------------
 

Tests/Utils/testChebInterpolator.cpp

@@ -104,7 +104,7 @@ int main(int, char **){
 	// approximative computation
 	const unsigned int ORDER = 10;
 	const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
-	typedef FChebInterpolator<ORDER> InterpolatorClass;
+	typedef FChebInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
 	InterpolatorClass S;
 
 	std::cout << "\nCompute interactions approximatively, interpolation order = " << ORDER << " ..." << std::endl;

Tests/Utils/testChebSxUCBSy.cpp

@@ -106,7 +106,7 @@ int main(int argc, char* argv[])
 	// approximative computation
 	const unsigned int ORDER = 10;
 	const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
-	typedef FChebInterpolator<ORDER> InterpolatorClass;
+	typedef FChebInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
 	InterpolatorClass S;
 	
 

UTests/utestChebyshevDirect.cpp

@@ -153,30 +153,30 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
 
 		// Print for information
 		Print("Potential diff is = ");
-		Print(potentialDiff.getL2Norm());
-		Print(potentialDiff.getInfNorm());
+		Print(potentialDiff.getRelativeL2Norm());
+		Print(potentialDiff.getRelativeInfNorm());
 		Print("Fx diff is = ");
-		Print(fx.getL2Norm());
-		Print(fx.getInfNorm());
+		Print(fx.getRelativeL2Norm());
+		Print(fx.getRelativeInfNorm());
 		Print("Fy diff is = ");
-		Print(fy.getL2Norm());
-		Print(fy.getInfNorm());
+		Print(fy.getRelativeL2Norm());
+		Print(fy.getRelativeInfNorm());
 		Print("Fz diff is = ");
-		Print(fz.getL2Norm());
-		Print(fz.getInfNorm());
+		Print(fz.getRelativeL2Norm());
+		Print(fz.getRelativeInfNorm());
 
 		// Assert
         const FReal MaximumDiffPotential = FReal(9e-5);
         const FReal MaximumDiffForces = FReal(9e-3);
 
-		uassert(potentialDiff.getL2Norm() < MaximumDiffPotential);
-		uassert(potentialDiff.getInfNorm() < MaximumDiffPotential);
-		uassert(fx.getL2Norm()  < MaximumDiffForces);
-		uassert(fx.getInfNorm() < MaximumDiffForces);
-		uassert(fy.getL2Norm()  < MaximumDiffForces);
-		uassert(fy.getInfNorm() < MaximumDiffForces);
-		uassert(fz.getL2Norm()  < MaximumDiffForces);