Introduced tensorial kernels for Uniform FMM (NB: Chebyshev UCB incompatible,...

Introduced tensorial kernels for Uniform FMM (NB: Chebyshev UCB incompatible, TODO Chebyshev Sym), tested for kernel R_{,ij} and its derivative, fixed multidimensionnal fft, + some cleanup.

Src/Kernels/Chebyshev/FChebCell.hpp

@@ -31,18 +31,18 @@
 * This class defines a cell used in the Chebyshev based FMM.
 * @param NVALS is the number of right hand side.
 */
-template <int ORDER, int NVALS = 1>
+template <int ORDER, int NMUL = 1, int NLOC = 1>
 class FChebCell : public FExtendMortonIndex, public FExtendCoordinate
 {
     static const int VectorSize = TensorTraits<ORDER>::nnodes * 2;
 
-    FReal multipole_exp[NVALS * VectorSize]; //< Multipole expansion
-    FReal     local_exp[NVALS * VectorSize]; //< Local expansion
+    FReal multipole_exp[NMUL * VectorSize]; //< Multipole expansion
+    FReal     local_exp[NLOC * VectorSize]; //< Local expansion
 	
 public:
     FChebCell(){
-        memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
-        memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
+        memset(multipole_exp, 0, sizeof(FReal) * NMUL * VectorSize);
+        memset(local_exp, 0, sizeof(FReal) * NLOC * VectorSize);
     }
 
 	~FChebCell() {}
@@ -72,26 +72,26 @@ public:
 
     /** Make it like the begining */
     void resetToInitialState(){
-        memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
-        memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
+        memset(multipole_exp, 0, sizeof(FReal) * NMUL * VectorSize);
+        memset(local_exp, 0, sizeof(FReal) * NLOC * VectorSize);
     }
 };
 
-template <int ORDER, int NVALS = 1>
-class FTypedChebCell : public FChebCell<ORDER,NVALS>, public FExtendCellType {
+template <int ORDER, int NMUL = 1, int NLOC = 1>
+class FTypedChebCell : public FChebCell<ORDER,NMUL,NLOC>, public FExtendCellType {
 public:
     template <class BufferWriterClass>
     void save(BufferWriterClass& buffer) const{
-        FChebCell<ORDER,NVALS>::save(buffer);
+        FChebCell<ORDER,NMUL,NLOC>::save(buffer);
         FExtendCellType::save(buffer);
     }
     template <class BufferReaderClass>
     void restore(BufferReaderClass& buffer){
-        FChebCell<ORDER,NVALS>::restore(buffer);
+        FChebCell<ORDER,NMUL,NLOC>::restore(buffer);
         FExtendCellType::restore(buffer);
     }
     void resetToInitialState(){
-        FChebCell<ORDER,NVALS>::resetToInitialState();
+        FChebCell<ORDER,NMUL,NLOC>::resetToInitialState();
         FExtendCellType::resetToInitialState();
     }
 };

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 "./FChebTensorialM2LHandler.hpp"
+
+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.
+ *
+ * PB: TODO 
+ * Erase this version since UCB and tensorial kernels are incompatible!
+ * Keep it until symmetric version is implemented.
+ *
+ * @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 {nK   = MatrixKernelClass::NK,
+        nRhs = MatrixKernelClass::NRHS,
+        nLhs = MatrixKernelClass::NLHS};
+
+	// private types
+	typedef FChebTensorialM2LHandler<ORDER,MatrixKernelClass> 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 FPoint& inBoxCenter,
+							const FReal inBoxWidth,
+							const FReal Epsilon)
+        : FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
+																																															inBoxCenter,
+																																															inBoxWidth),
+			M2LHandler(new M2LHandlerClass(Epsilon))
+	{
+		// read precomputed compressed m2l operators from binary file
+		//M2LHandler->ReadFromBinaryFileAndSet();
+		M2LHandler->ComputeAndCompressAndSet();
+	}
+
+
+	void P2M(CellClass* const LeafCell,
+					 const ContainerClass* const SourceParticles)
+	{
+        const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
+//    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+        for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
+
+            // 1) apply Sy
+            AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
+                                                      LeafCell->getMultipole(idxRhs), SourceParticles);
+            // 2) apply B
+            M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes, idxRhs);
+        }
+//  }//NVALS
+	}
+
+
+	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){
+
+            // 1) apply Sy
+            FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxRhs));
+            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
+                if (ChildCells[ChildIndex]){
+                    AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
+                                                              ParentCell->getMultipole(idxRhs));
+                }
+            }
+            // 2) apply B
+            M2LHandler->applyB(ParentCell->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes, idxRhs);
+        }
+//  }//NVALS
+	}
+
+
+//	void M2L(CellClass* const FRestrict TargetCell,
+//					 const CellClass* SourceCells[343],
+//					 const int NumSourceCells,
+//					 const int TreeLevel) const
+//	{
+//		const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
+//		const FTreeCoordinate& cx = TargetCell->getCoordinate();
+//		for (int idx=0; idx<NumSourceCells; ++idx) {
+//			const FTreeCoordinate& cy = SourceCells[idx]->getCoordinate();
+//			const int transfer[3] = {cy.getX()-cx.getX(),
+//															 cy.getY()-cx.getY(),
+//															 cy.getZ()-cx.getZ()};
+//			M2LHandler->applyC(transfer, CellWidth,
+//												SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
+//												TargetCell->getLocal() + AbstractBaseClass::nnodes);
+//		}
+//	}
+
+	void M2L(CellClass* const FRestrict TargetCell,
+					 const CellClass* SourceCells[343],
+                     const int /*NumSourceCells*/,
+					 const int TreeLevel)
+	{
+//    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+    for (unsigned int idxLhs=0; idxLhs < nLhs; ++idxLhs)
+      for (unsigned int idxRhs=0; idxRhs < nRhs; ++idxRhs){
+        unsigned int idxK = idxLhs*nRhs + idxRhs;
+
+        FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxLhs) + AbstractBaseClass::nnodes;
+        const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
+        for (int idx=0; idx<343; ++idx){
+          if (SourceCells[idx]){
+            M2LHandler->applyC(idx, CellWidth, SourceCells[idx]->getMultipole(idxRhs) + AbstractBaseClass::nnodes,
+                               CompressedLocalExpansion, idxK);
+          }
+        }
+      }
+//  }//NVALS
+	}
+
+//	void M2L(CellClass* const FRestrict TargetCell,
+//					 const CellClass* SourceCells[343],
+//					 const int NumSourceCells,
+//					 const int TreeLevel) const
+//	{
+//		const unsigned int rank = M2LHandler.getRank();
+//		FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
+//		const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
+//		for (int idx=0; idx<343; ++idx)
+//			if (SourceCells[idx])
+//				FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idx]->getMultipole())+AbstractBaseClass::nnodes,
+//										MultipoleExpansion+idx*rank);
+//		
+//		M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
+//	}
+
+
+	void L2L(const CellClass* const FRestrict ParentCell,
+					 CellClass* FRestrict *const FRestrict ChildCells,
+                     const int /*TreeLevel*/)
+	{
+//    for(int idxV = 0 ; idxV < NVALS ; ++idxV){
+    for(int idxRhs = 0 ; idxRhs < nLhs ; ++idxRhs){
+
+            // 1) apply U
+            M2LHandler->applyU(ParentCell->getLocal(idxRhs) + AbstractBaseClass::nnodes,
+                               const_cast<CellClass*>(ParentCell)->getLocal(idxRhs), idxRhs);
+            // 2) apply Sx
+            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
+                if (ChildCells[ChildIndex]){
+                    AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
+                }
+            }
+        }
+//  }//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 idxRhs = 0 ; idxRhs < nLhs ; ++idxRhs){
+
+            // 1) apply U
+          M2LHandler->applyU(LeafCell->getLocal(idxRhs) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxRhs), idxRhs);
+
+            //// 2.a) apply Sx
+            //AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
+            //																					AbstractBaseClass::BoxWidthLeaf,
+            //																					LeafCell->getLocal(),
+            //																					TargetParticles);
+            //// 2.b) apply Px (grad Sx)
+            //AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
+            //																									AbstractBaseClass::BoxWidthLeaf,
+            //																									LeafCell->getLocal(),
+            //																									TargetParticles);
+
+            // 2.c) apply Sx and Px (grad Sx)
+            AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
+                                                            LeafCell->getLocal(idxRhs), TargetParticles);
+        }
+//  }//NVALS
+	}
+
+    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

@@ -16,16 +16,20 @@
 #ifndef FINTERPMATRIXKERNEL_HPP
 #define FINTERPMATRIXKERNEL_HPP
 
+#include <stdexcept>
+
 #include "../../Utils/FPoint.hpp"
 #include "../../Utils/FNoCopyable.hpp"
 #include "../../Utils/FMath.hpp"
 #include "../../Utils/FBlas.hpp"
 
-
 // extendable
 enum KERNEL_FUNCTION_IDENTIFIER {ONE_OVER_R,
                                  ONE_OVER_R_SQUARED,
-                                 LENNARD_JONES_POTENTIAL};
+                                 LENNARD_JONES_POTENTIAL,
+                                 ID_OVER_R,
+                                 R_IJ,
+                                 R_IJK};
 
 // probably not extedable :)
 enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
@@ -52,9 +56,35 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
 {
   static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R;
+  static const  unsigned int Dim = 1; //PB: dimension of kernel
+  const unsigned int indexTab[2]={0,
+                                  0};
+  static const unsigned int NK = 1;
+  static const unsigned int NRHS = 1;
+  static const unsigned int NLHS = 1;
+
+  // PB: figure out if applyTab is already of size NRHS*NLHS ?
+  const unsigned int applyTab[1]={0}; // get position in sym tensor
+
+  unsigned int _i,_j;
 
   FInterpMatrixKernelR() {}
 
+  // updates indices _i and _j for current position d in reduced storage 
+  void updateIndex(const unsigned int d) //const
+  {_i=indexTab[d]; _j=indexTab[d+Dim];}
+
+  // returns indices i and j for a given position in full storage
+  int getIndexLhs(const unsigned int d) const
+  {return indexTab[applyTab[d]];}
+
+  int getIndexRhs(const unsigned int d) const
+  {return indexTab[applyTab[d]+Dim];}
+
+  // returns position in reduced storage (from n= position in full ?or? for i, j)
+  int getPosition(const unsigned int n) const
+  {return applyTab[n];} 
+
   FReal evaluate(const FPoint& x, const FPoint& y) const
   {
     const FPoint xy(x-y);
@@ -82,7 +112,15 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
 {
   static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
-
+  static const  unsigned int Dim = 1; //PB: dimension of kernel
+  const unsigned int indexTab[2]={0,
+                                  0};
+  static const unsigned int NK = 1;
+  static const unsigned int NRHS = 1;
+  static const unsigned int NLHS = 1;
+
+  // PB: figure out if applyTab is already of size NRHS*NLHS ?
+  const unsigned int applyTab[1]={0}; // get position in sym tensor
   FInterpMatrixKernelRR() {}
 
   FReal evaluate(const FPoint& x, const FPoint& y) const
@@ -112,6 +150,7 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
 {
   static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
   static const KERNEL_FUNCTION_IDENTIFIER Identifier = LENNARD_JONES_POTENTIAL;
+  static const  unsigned int Dim = 1; //PB: dimension of kernel
 
   FInterpMatrixKernelLJ() {}
 
@@ -137,9 +176,236 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
     // return 1 because non homogeneous kernel functions cannot be scaled!!!
     return FReal(1.);
   }
+
 };
 
 
+/// Test Tensorial kernel 1/R*Id_3
+struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
+{
+  static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
+  static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R;
+  static const unsigned int Dim = 6; //PB: dimension of kernel
+  const unsigned int indexTab[12]={0,0,0,1,1,2,
+                                   0,1,2,1,2,2};
+  static const unsigned int NK = 9;
+  static const unsigned int NRHS = 3;
+  static const unsigned int NLHS = 3;
+
+  // PB: figure out if applyTab is already of size NRHS*NLHS ?
+  const unsigned int applyTab[9]={0,1,2,1,3,4,2,4,5}; // get position in sym tensor
+
+  unsigned int _i,_j;
+
+  FInterpMatrixKernel_IOR() {}
+
+  // updates indices _i and _j for current position d in reduced storage 
+  void updateIndex(const unsigned int d) //const
+  {_i=indexTab[d]; _j=indexTab[d+Dim];}
+
+  // returns indices i and j for a given position in full storage
+  int getIndexLhs(const unsigned int d) const
+  {return indexTab[applyTab[d]];}
+
+  int getIndexRhs(const unsigned int d) const
+  {return indexTab[applyTab[d]+Dim];}
+
+  // returns position in reduced storage (from n= position in full ?or? for i, j)
+  int getPosition(const unsigned int n) const
+  {return applyTab[n];} 
+
+  FReal evaluate(const FPoint& x, const FPoint& y) const
+  {
+    const FPoint xy(x-y);
+
+    if(_i==_j)
+      return FReal(1.)/xy.norm();
+    else
+      return FReal(0.);
+
+  }
+
+  FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
+  {
+    const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
+    return getScaleFactor(CellWidth);
+  }
+
+  FReal getScaleFactor(const FReal CellWidth) const
+  {
+        return FReal(2.) / CellWidth;
+  }
+
+};
+
+/// R_{,ij}
+struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
+{
+  static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
+  static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJ;
+  static const unsigned int Dim = 6; //PB: dimension of kernel
+  const unsigned int indexTab[12]={0,0,0,1,1,2,
+                                   0,1,2,1,2,2};
+  static const unsigned int NK = 9;
+  static const unsigned int NRHS = 3;
+  static const unsigned int NLHS = 3;
+
+  // PB: figure out if applyTab is already of size NRHS*NLHS ?
+  const unsigned int applyTab[9]={0,1,2,1,3,4,2,4,5}; // get position in sym tensor
+
+  unsigned int _i,_j;
+
+  FInterpMatrixKernel_R_IJ() {}
+
+  // updates indices _i and _j for current position d in reduced storage 
+  void updateIndex(const unsigned int d) //const
+  {_i=indexTab[d]; _j=indexTab[d+Dim];}
+
+  // returns indices i and j for a given position in full storage
+  int getIndexLhs(const unsigned int d) const
+  {return indexTab[applyTab[d]];}
+
+  int getIndexRhs(const unsigned int d) const
+  {return indexTab[applyTab[d]+Dim];}
+
+  // returns position in reduced storage (from n= position in full ?or? for i, j)
+  int getPosition(const unsigned int n) const
+  {return applyTab[n];} 
+
+  FReal evaluate(const FPoint& x, const FPoint& y) const
+  {
+    const FPoint xy(x-y);
+    const FReal one_over_r = FReal(1.)/xy.norm();
+    const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
+    double ri,rj;
+
+    if(_i==0) ri=xy.getX();
+    else if(_i==1) ri=xy.getY();
+    else if(_i==2) ri=xy.getZ();
+    else throw std::runtime_error("Update i!");
+
+    if(_j==0) rj=xy.getX();
+    else if(_j==1) rj=xy.getY();
+    else if(_j==2) rj=xy.getZ();
+    else throw std::runtime_error("Update j!");
+
+    if(_i==_j)
+      return one_over_r - ri * ri * one_over_r3;
+    else
+      return - ri * rj * one_over_r3;
+
+  }
+
+  FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
+  {
+    const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
+    return getScaleFactor(CellWidth);
+  }
+
+  // R_{,ij} is homogeneous to [L]/[L]^{-2}=[L]^{-1}
+  // => same scale factor as ONE_OVER_R
+  FReal getScaleFactor(const FReal CellWidth) const
+  {
+        return FReal(2.) / CellWidth;
+  }
+
+};
+
+
+/// R_{,ijk}
+struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
+{
+  static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
+  static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJK;
+  static const  unsigned int Dim = 10; //PB: dimension of kernel
+  const unsigned int indexTab[30]={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 NK = 27;
+  static const unsigned int NRHS = 3;
+  static const unsigned int NLHS = 3;
+
+  // PB: figure out if applyTab is already of size NRHS*NLHS ?
+  const unsigned int applyTab[27]={0,3,6,3,1,9,6,9,2,
+                                   3,1,9,1,4,7,9,7,5,
+                                   6,9,2,9,7,5,2,5,8}; // get position in sym tensor
+
+
+  unsigned int _i,_j,_k;
+
+  FInterpMatrixKernel_R_IJK() {}
+
+  void updateIndex(unsigned int d) //const
+  {
+    _i=indexTab[d]; _j=indexTab[d+Dim]; _k=indexTab[d+2*Dim];
+  }
+
+  // returns indices i and j for a given position in full storage
+  int getIndexLhs(const unsigned int d) const
+  {return indexTab[applyTab[d]];}
+
+  int getIndexRhs(const unsigned int d) const
+  {return indexTab[applyTab[d]+Dim];}
+
+  // returns position in reduced storage (from n= position in full ?or? for i, j)
+  int getPosition(const unsigned int n) const
+  {return applyTab[n];} 
+
+  FReal evaluate(const FPoint& x, const FPoint& y) const
+  {
+    const FPoint xy(x-y);
+    const FReal one_over_r = FReal(1.)/xy.norm();
+    const FReal one_over_r2 = one_over_r*one_over_r;
+    const FReal one_over_r3 = one_over_r2*one_over_r;
+    double ri,rj,rk;
+
+    if(_i==0) ri=xy.getX();
+    else if(_i==1) ri=xy.getY();
+    else if(_i==2) ri=xy.getZ();
+    else throw std::runtime_error("Update i!");
+
+    if(_j==0) rj=xy.getX();
+    else if(_j==1) rj=xy.getY();
+    else if(_j==2) rj=xy.getZ();
+    else throw std