Fixed getInteractionNeighbors to work with new separation criteria; Propagated...

Fixed getInteractionNeighbors to work with new separation criteria; Propagated separation criterion into Uniform kernel and m2lhandler (for precomputation); Provided test with one over r2 with corewidth (pb with derivative).

Src/Containers/FOctree.hpp

@@ -832,7 +832,7 @@ public:
                     if(!FMath::Between(parentCell.getZ() + idxZ,0,boxLimite)) continue;
 
                     // if we are not on the current cell
-                    if( idxX || idxY || idxZ ){
+                    if( neighSeparation<1 || 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

Src/Core/FFmmAlgorithm.hpp

@@ -47,7 +47,7 @@ class FFmmAlgorithm :  public FAbstractAlgorithm, public FAlgorithmTimers {
 
     const int OctreeHeight;        ///< The height of the given tree.
 
-    const int leafLevelSeperationCriteria;
+    const int leafLevelSeparationCriteria;
 public:
     /** Class constructor
      * 
@@ -57,8 +57,8 @@ public:
      *
      * \except An exception is thrown if one of the arguments is NULL.
      */
-    FFmmAlgorithm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeperationCriteria = 1)
-        : tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()), leafLevelSeperationCriteria(inLeafLevelSeperationCriteria) {
+    FFmmAlgorithm(OctreeClass* const inTree, KernelClass* const inKernels, const int inLeafLevelSeparationCriteria = 1)
+        : tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()), leafLevelSeparationCriteria(inLeafLevelSeparationCriteria) {
 
         FAssertLF(tree, "tree cannot be null");
         FAssertLF(kernels, "kernels cannot be null");
@@ -191,12 +191,11 @@ protected:
 
         const CellClass* neighbors[343];
 
-
         // for each levels
         for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
             FLOG(FTic counterTimeLevel);
 
-            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeperationCriteria);
+            const int separationCriteria = (idxLevel != FAbstractAlgorithm::lowerWorkingLevel-1 ? 1 : leafLevelSeparationCriteria);
 
             // for each cells
             do{

Src/Kernels/Interpolation/FInterpMatrixKernel.hpp

@@ -433,6 +433,107 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel<FReal>
     }
 };
 
+
+/// One over (a+r^2)
+template <class FReal>
+struct FInterpMatrixKernelAPLUSRR : FInterpAbstractMatrixKernel<FReal>
+{
+    static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
+    static const unsigned int NCMP = 1; //< number of components
+    static const unsigned int NPV  = 1; //< dim of physical values
+    static const unsigned int NPOT = 1; //< dim of potentials
+    static const unsigned int NRHS = 1; //< dim of mult exp
+    static const unsigned int NLHS = 1; //< dim of loc exp
+
+    const FReal CoreWidth;
+
+    FInterpMatrixKernelAPLUSRR(const FReal inCoreWidth = .25)
+    : CoreWidth(inCoreWidth)
+     {}
+
+    // copy ctor
+    FInterpMatrixKernelAPLUSRR(const FInterpMatrixKernelAPLUSRR& other)
+    : CoreWidth(other.CoreWidth)
+    {}
+
+  static const char* getID() { return "ONE_OVER_A_PLUS_RR"; }
+
+    // returns position in reduced storage
+    int getPosition(const unsigned int) const
+    {return 0;}
+
+    // evaluate interaction
+    template <class ValueClass>
+    ValueClass evaluate(const ValueClass& x1, const ValueClass& y1, const ValueClass& z1,
+                   const ValueClass& x2, const ValueClass& y2, const ValueClass& z2) const
+    {
+        const ValueClass diffx = (x1-x2);
+        const ValueClass diffy = (y1-y2);
+        const ValueClass diffz = (z1-z2);
+        return FMath::One<ValueClass>() / (FMath::ConvertTo<ValueClass,FReal>(CoreWidth) + // WHY FReal??
+                                           diffx*diffx +
+                                           diffy*diffy +
+                                           diffz*diffz);
+    }
+
+    // evaluate interaction (blockwise)
+    template <class ValueClass>
+    void evaluateBlock(const ValueClass& x1, const ValueClass& y1, const ValueClass& z1,
+                       const ValueClass& x2, const ValueClass& y2, const ValueClass& z2, ValueClass* block) const
+    {
+        block[0]=this->evaluate(x1,y1,z1,x2,y2,z2);
+    }
+
+    // evaluate interaction and derivative (blockwise)
+    template <class ValueClass>
+    void evaluateBlockAndDerivative(const ValueClass& x1, const ValueClass& y1, const ValueClass& z1,
+                                    const ValueClass& x2, const ValueClass& y2, const ValueClass& z2,
+                                    ValueClass block[1], ValueClass blockDerivative[3]) const
+    {
+        const ValueClass diffx = (x1-x2);
+        const ValueClass diffy = (y1-y2);
+        const ValueClass diffz = (z1-z2);
+        const ValueClass r2 = (diffx*diffx +
+                               diffy*diffy +
+                               diffz*diffz);
+        const ValueClass one_over_a_plus_r2 = FMath::One<ValueClass>() / (FMath::ConvertTo<ValueClass,FReal>(CoreWidth)+r2);
+        const ValueClass one_over_a_plus_r2_squared = one_over_a_plus_r2*one_over_a_plus_r2;
+
+        block[0] = one_over_a_plus_r2;
+
+        // TODO Fix derivative
+        const ValueClass coef = FMath::ConvertTo<ValueClass,FReal>(-2.) * one_over_a_plus_r2_squared;
+        blockDerivative[0] = coef * diffx;
+        blockDerivative[1] = coef * diffy;
+        blockDerivative[2] = coef * diffz;
+
+    }
+
+    FReal getScaleFactor(const FReal, const int) const
+    {
+        // return 1 because non homogeneous kernel functions cannot be scaled!!!
+        return FReal(1.0);
+    }
+
+    FReal getScaleFactor(const FReal) const
+    {
+        // return 1 because non homogeneous kernel functions cannot be scaled!!!
+        return FReal(1.0);    
+    }
+
+    FReal evaluate(const FPoint<FReal>& p1, const FPoint<FReal>& p2) const{
+        return evaluate<FReal>(p1.getX(), p1.getY(), p1.getZ(), p2.getX(), p2.getY(), p2.getZ());
+    }
+    void evaluateBlock(const FPoint<FReal>& p1, const FPoint<FReal>& p2, FReal* block) const{
+        evaluateBlock<FReal>(p1.getX(), p1.getY(), p1.getZ(), p2.getX(), p2.getY(), p2.getZ(), block);
+    }
+    void evaluateBlockAndDerivative(const FPoint<FReal>& p1, const FPoint<FReal>& p2,
+                                    FReal block[1], FReal blockDerivative[3]) const {
+        evaluateBlockAndDerivative<FReal>(p1.getX(), p1.getY(), p1.getZ(), p2.getX(), p2.getY(), p2.getZ(), block, blockDerivative);
+    }
+};
+
+
 ////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////
 //

Src/Kernels/Uniform/FUnifKernel.hpp

@@ -60,6 +60,8 @@ class FUnifKernel
     /// Needed for M2L operator
     const M2LHandlerClass M2LHandler;
 
+    /// Leaf level separation criterion
+    const int LeafLevelSeparationCriterion;
 
 public:
     /**
@@ -70,12 +72,15 @@ public:
     FUnifKernel(const int inTreeHeight,
                 const FReal inBoxWidth,
                 const FPoint<FReal>& inBoxCenter,
-                const MatrixKernelClass *const inMatrixKernel)
+                const MatrixKernelClass *const inMatrixKernel,
+                const int inLeafLevelSeparationCriterion)
     : FAbstractUnifKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
       MatrixKernel(inMatrixKernel),
       M2LHandler(MatrixKernel,
                  inTreeHeight,
-                 inBoxWidth) 
+                 inBoxWidth,
+                 inLeafLevelSeparationCriterion) ,
+      LeafLevelSeparationCriterion(inLeafLevelSeparationCriterion)
     { }
 
 
@@ -186,7 +191,13 @@ public:
              ContainerClass* const NeighborSourceParticles[27],
              const int /* size */)
     {
-        DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
+        // Standard FMM separation criterion, i.e. max 27 neighbor clusters per leaf
+        if(LeafLevelSeparationCriterion==1) 
+            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
+        // Nearfield interactions are only computed within the target leaf
+        if(LeafLevelSeparationCriterion==0) 
+            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(TargetParticles,NeighborSourceParticles,0,MatrixKernel);
+        // If criterion equals -1 then no P2P need to be performed.
     }
 
 
@@ -194,7 +205,13 @@ public:
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/)
     {
-        DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
+        // Standard FMM separation criterion, i.e. max 27 neighbor clusters per leaf
+        if(LeafLevelSeparationCriterion==1) 
+            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
+        // Nearfield interactions are only computed within the target leaf
+        if(LeafLevelSeparationCriterion==0) 
+            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,0,MatrixKernel);
+        // If criterion equals -1 then no P2P need to be performed.
     }
 
 };

Src/Kernels/Uniform/FUnifM2LHandler.hpp

@@ -38,14 +38,14 @@
 /*!  Precomputation of the 316 interactions by evaluation of the matrix kernel on the uniform grid and transformation into Fourier space. 
 PB: Compute() does not belong to the M2LHandler like it does in the Chebyshev kernel. This allows much nicer specialization of the M2LHandler class with respect to the homogeneity of the kernel of interaction like in the ChebyshevSym kernel.*/
 template < class FReal,int ORDER, typename MatrixKernelClass>
-static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, FComplex<FReal>* &FC)
+static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth, FComplex<FReal>* &FC, const int SeparationCriterion = 1)
 {
     // allocate memory and store compressed M2L operators
     if (FC) throw std::runtime_error("M2L operators are already set");
     // dimensions of operators
     const unsigned int order = ORDER;
     const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
-    const unsigned int ninteractions = 316;
+    const unsigned int ninteractions = 316+26*(SeparationCriterion<1 ? 1 : 0) + 1*(SeparationCriterion<0 ? 1 : 0);
     typedef FUnifTensor<FReal,ORDER> TensorType;
 
     // interpolation points of source (Y) and target (X) cell
@@ -53,7 +53,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
     // set roots of target cell (X)
     TensorType::setRoots(FPoint<FReal>(0.,0.,0.), CellWidth, X);
 
-    // allocate memory and compute 316 m2l operators
+    // allocate memory and compute 316 m2l operators (342 if separation equals 0, 343 if separation equals -1)
     FReal    *_C;
     FComplex<FReal> *_FC;
 
@@ -78,7 +78,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
     for (int i=-3; i<=3; ++i) {
         for (int j=-3; j<=3; ++j) {
             for (int k=-3; k<=3; ++k) {
-                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
+                if (abs(i)>SeparationCriterion || abs(j)>SeparationCriterion || abs(k)>SeparationCriterion) {
                     // set roots of source cell (Y)
                     const FPoint<FReal> cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
                     FUnifTensor<FReal,order>::setRoots(cy, CellWidth, Y);
@@ -110,7 +110,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
         }
     }
     if (counter != ninteractions)
-        throw std::runtime_error("Number of interactions must correspond to 316");
+        throw std::runtime_error("Number of interactions must correspond to " + std::to_string(ninteractions));
 
     // Free _C
     delete [] _C;
@@ -126,7 +126,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
         for (int j=-3; j<=3; ++j)
             for (int k=-3; k<=3; ++k) {
                 const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
-                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
+                if (abs(i)>SeparationCriterion || abs(j)>SeparationCriterion || abs(k)>SeparationCriterion) {
                     FBlas::c_copy(opt_rc, reinterpret_cast<FReal*>(_FC + counter*rc), 
                                   reinterpret_cast<FReal*>(FC + idx*opt_rc));
                     counter++;
@@ -136,7 +136,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
       }
 
     if (counter != ninteractions)
-        throw std::runtime_error("Number of interactions must correspond to 316");
+        throw std::runtime_error("Number of interactions must correspond to " + std::to_string(ninteractions));
     delete [] _FC;      
 }
 
@@ -182,6 +182,8 @@ class FUnifM2LHandler<FReal, ORDER,HOMOGENEOUS>
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
+    /// Leaf level separation criterion
+    const int LeafLevelSeparationCriterion;
 
     static const std::string getFileName()
     {
@@ -195,8 +197,8 @@ class FUnifM2LHandler<FReal, ORDER,HOMOGENEOUS>
     
 public:
     template <typename MatrixKernelClass>
-    FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal)
-        : FC(nullptr), Dft(), opt_rc(rc/2+1)
+    FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const int inLeafLevelSeparationCriterion)
+        : FC(nullptr), Dft(), opt_rc(rc/2+1), LeafLevelSeparationCriterion(inLeafLevelSeparationCriterion)
     {    
         // init DFT
         const int steps[dimfft] = {rc};
@@ -239,7 +241,7 @@ public:
         // Compute matrix of interactions
         const FReal ReferenceCellWidth = FReal(2.);
         FComplex<FReal>* pFC = NULL;
-        Compute<FReal,order>(MatrixKernel,ReferenceCellWidth,pFC);
+        Compute<FReal,order>(MatrixKernel,ReferenceCellWidth,pFC,LeafLevelSeparationCriterion);
         FC.assign(pFC);
 
         // Compute memory usage
@@ -250,7 +252,10 @@ public:
         std::cout << "Compute and set M2L operators ("<< long(sizeM2L/**1e-6*/) <<" B) in "
                   << time.tacAndElapsed() << "sec."   << std::endl;
     }
-        
+
+    unsigned long long getMemory() const {
+        return 343*opt_rc*sizeof(FComplex<FReal>);
+    }        
 
     /**
     * Expands potentials \f$x+=IDFT(X)\f$ of a target cell. This operation can be
@@ -341,6 +346,8 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
+    /// Leaf level separation criterion
+    const int LeafLevelSeparationCriterion;
 
     static const std::string getFileName()
     {
@@ -354,10 +361,10 @@ class FUnifM2LHandler<FReal,ORDER,NON_HOMOGENEOUS>
     
 public:
     template <typename MatrixKernelClass>
-    FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth)
+    FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const int inLeafLevelSeparationCriterion)
         : TreeHeight(inTreeHeight),
           RootCellWidth(inRootCellWidth),
-          Dft(), opt_rc(rc/2+1)
+          Dft(), opt_rc(rc/2+1), LeafLevelSeparationCriterion(inLeafLevelSeparationCriterion)
     {
         // init DFT
         const int steps[dimfft] = {rc};
@@ -414,9 +421,12 @@ public:
         CellWidth /= FReal(2.);                      // at level 2
         for (unsigned int l=2; l<TreeHeight; ++l) {
 
+            // Determine separation criteria wrt level
+            const int SeparationCriterion = (l != TreeHeight-1 ? 1 : LeafLevelSeparationCriterion);
+
             // check if already set
             if (FC[l]) throw std::runtime_error("M2L operator already set");
-            Compute<FReal,order>(MatrixKernel,CellWidth,FC[l]);
+            Compute<FReal,order>(MatrixKernel,CellWidth,FC[l],SeparationCriterion);
             CellWidth /= FReal(2.);                    // at level l+1 
 
         }
@@ -429,6 +439,9 @@ public:
                   << time.tacAndElapsed() << "sec."   << std::endl;
     }
 
+    unsigned long long getMemory() const {
+        return (TreeHeight-2)*343*opt_rc*sizeof(FComplex<FReal>);
+    }   
 
     /**
     * Expands potentials \f$x+=IDFT(X)\f$ of a target cell. This operation can be

Src/Kernels/Uniform/FUnifTensorialM2LHandler.hpp

@@ -42,12 +42,13 @@ template < class FReal, int ORDER, class MatrixKernelClass>
 static void Compute(const MatrixKernelClass *const MatrixKernel, 
                     const FReal CellWidth, 
                     const FReal CellWidthExtension, 
-                    FComplex<FReal>** &FC)
+                    FComplex<FReal>** &FC, 
+                    const int SeparationCriterion = 1)
 {
     // dimensions of operators
     const unsigned int order = ORDER;
     const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
-    const unsigned int ninteractions = 316;
+    const unsigned int ninteractions = 316+26*(SeparationCriterion<1 ? 1 : 0) + 1*(SeparationCriterion<0 ? 1 : 0);
     const unsigned int ncmp = MatrixKernelClass::NCMP;
 
     // utils
@@ -90,11 +91,14 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
     unsigned int perm[rc];
     TensorType::setStoragePermutation(perm);
 
+    // Allocate intermediate interaction block
+    FReal* block = new FReal[ncmp]; 
+
     unsigned int counter = 0;
     for (int i=-3; i<=3; ++i) {
         for (int j=-3; j<=3; ++j) {
             for (int k=-3; k<=3; ++k) {
-                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
+                if (abs(i)>SeparationCriterion || abs(j)>SeparationCriterion || abs(k)>SeparationCriterion) {
                     // set roots of source cell (Y)
                     const FPoint<FReal> cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
                     FUnifTensor<FReal,order>::setRoots(cy, ExtendedCellWidth, Y);
@@ -104,8 +108,6 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
                         for(unsigned int m=0; m<2*order-1; ++m)
                             for(unsigned int n=0; n<2*order-1; ++n){   
                                 // Compute current M2L interaction (block matrix)
-                                FReal* block; 
-                                block = new FReal[ncmp]; 
                                 MatrixKernel->evaluateBlock(X[node_ids_pairs[ido][0]], 
                                                             Y[node_ids_pairs[ido][1]],
                                                             block);
@@ -133,6 +135,9 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
     if (counter != ninteractions)
         throw std::runtime_error("Number of interactions must correspond to 316");
 
+    // Free block
+    delete [] block;
+
     // Free _C
     for (unsigned int d=0; d<ncmp; ++d) 
     delete [] _C[d];
@@ -149,7 +154,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
         for (int j=-3; j<=3; ++j)
             for (int k=-3; k<=3; ++k) {
                 const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
-                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
+                if (abs(i)>SeparationCriterion || abs(j)>SeparationCriterion || abs(k)>SeparationCriterion) {
                     for (unsigned int d=0; d<ncmp; ++d) 
                         FBlas::c_copy(opt_rc, reinterpret_cast<FReal*>(_FC[d] + counter*rc), 
                                       reinterpret_cast<FReal*>(FC[d] + idx*opt_rc));
@@ -207,6 +212,8 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
+    /// Leaf level separation criterion
+    const int LeafLevelSeparationCriterion;
 
     static const std::string getFileName()
     {
@@ -219,9 +226,9 @@ class FUnifTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS>
 
     
 public:
-    FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal inCellWidthExtension)
+    FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal inCellWidthExtension, const int inLeafLevelSeparationCriterion)
         :  CellWidthExtension(inCellWidthExtension), 
-           Dft(), opt_rc(rc/2+1)
+           Dft(), opt_rc(rc/2+1), LeafLevelSeparationCriterion(inLeafLevelSeparationCriterion)
     {
         // init DFT
         const int steps[dimfft] = {rc};
@@ -245,7 +252,7 @@ public:
   FUnifTensorialM2LHandler(const FUnifTensorialM2LHandler& other)
     : FC(other.FC),
       CellWidthExtension(other.CellWidthExtension),
-      Dft(), opt_rc(other.opt_rc)
+      Dft(), opt_rc(other.opt_rc), LeafLevelSeparationCriterion(other.LeafLevelSeparationCriterion)
     {    
         // init DFT
         const int steps[dimfft] = {rc};
@@ -278,7 +285,7 @@ public:
         // but it NEEDS to match the numerator of the scale factor in matrix kernel!
         // Therefore box width extension is not yet supported for homog kernels
         const FReal ReferenceCellWidth = FReal(2.);
-        Compute<FReal,order>(MatrixKernel,ReferenceCellWidth, 0., FC);
+        Compute<FReal,order>(MatrixKernel,ReferenceCellWidth, 0., FC, LeafLevelSeparationCriterion);
         
         // Compute memory usage
         unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
@@ -288,6 +295,10 @@ public:
                   << time.tacAndElapsed() << "sec."   << std::endl;
     }
 
+    unsigned long long getMemory() const {
+        return 343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+    }        
+
     /**
      * Expands potentials \f$x+=IDFT(X)\f$ of a target cell. This operation can be
      * seen as part of the L2L operation.
@@ -396,6 +407,8 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
     DftClass Dft;
     const unsigned int opt_rc; // specific to real valued kernel
 
+    /// Leaf level separation criterion
+    const int LeafLevelSeparationCriterion;
 
     static const std::string getFileName()
     {
@@ -408,11 +421,11 @@ class FUnifTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
 
     
 public:
-    FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal inCellWidthExtension)
+    FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal inCellWidthExtension, const int inLeafLevelSeparationCriterion)
         : TreeHeight(inTreeHeight),
           RootCellWidth(inRootCellWidth), 
           CellWidthExtension(inCellWidthExtension),
-          Dft(), opt_rc(rc/2+1)
+          Dft(), opt_rc(rc/2+1), LeafLevelSeparationCriterion(inLeafLevelSeparationCriterion)
     {
         // init DFT
         const int steps[dimfft] = {rc};
@@ -441,7 +454,7 @@ public:
       TreeHeight(other.TreeHeight),
       RootCellWidth(other.RootCellWidth), 
       CellWidthExtension(other.CellWidthExtension),
-      Dft(), opt_rc(other.opt_rc)
+      Dft(), opt_rc(other.opt_rc), LeafLevelSeparationCriterion(other.LeafLevelSeparationCriterion)
     {    
         // init DFT
         const int steps[dimfft] = {rc};
@@ -471,10 +484,12 @@ public:
         FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
         CellWidth /= FReal(2.);                      // at level 2
         for (unsigned int l=2; l<TreeHeight; ++l) {
+            // Determine separation criteria wrt level
+            const int SeparationCriterion = (l != TreeHeight-1 ? 1 : LeafLevelSeparationCriterion);
             // check if already set
             for (unsigned int d=0; d<ncmp; ++d)
                 if (FC[l][d]) throw std::runtime_error("M2L operator already set");
-            Compute<FReal,order>(MatrixKernel,CellWidth,CellWidthExtension,FC[l]);
+            Compute<FReal,order>(MatrixKernel,CellWidth,CellWidthExtension,FC[l],SeparationCriterion);
             CellWidth /= FReal(2.);                    // at level l+1 
         }
 
@@ -486,6 +501,10 @@ public:
                   << time.tacAndElapsed() << "sec."   << std::endl;
     }
 
+    unsigned long long getMemory() const {
+        return (TreeHeight-2)*343*ncmp*opt_rc*sizeof(FComplex<FReal>);
+    }        
+
     /**
      * Expands potentials \f$x+=IDFT(X)\f$ of a target cell. This operation can be
      * seen as part of the L2L operation.

Src/Utils/FParameterNames.hpp

@@ -182,6 +182,11 @@ static const FParameterNames PhysicalValue = {
      "The physical value of the particles."
 };
 
+static const FParameterNames SeparationCriterion = {
+    {"-sep", "--separation-criterion"} ,
+     "Specify number of clusters separing 2 well-separated clusters."
+};
+
 /** To print a list of parameters */
 inline void PrintUsedOptions(const std::vector<FParameterNames>& options){
     std::cout << ">> Here is the list of the parameters you can pass to this application :\n";

Tests/Kernels/testSmoothUnifAlgorithm.cpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger 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".
+// ===================================================================================
+
+
+/**
+ *@author Pierre Blanchard
+ *
+ * **/
+// ==== CMAKE =====
+// @FUSE_FFT
+// ================
+// Keep in private GIT
+// @SCALFMM_PRIVATE
+
+
+#include <iostream>
+
+#include <cstdio>
+#include <cstdlib>
+
+#include "Files/FFmaGenericLoader.hpp"
+
+
+#include "Kernels/Uniform/FUnifCell.hpp"
+#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
+#include "Kernels/Uniform/FUnifKernel.hpp"
+
+#include "Components/FSimpleLeaf.hpp"
+#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
+
+#include "Utils/FParameters.hpp"
+#include "Utils/FMemUtils.hpp"
+
+#include "Containers/FOctree.hpp"
+#include "Containers/FVector.hpp"
+
+#include "Core/FFmmAlgorithm.hpp"
+#include "Core/FFmmAlgorithmThread.hpp"
+
+#include "../../Src/Utils/FParameterNames.hpp"
+
+/**
+ * This program runs the FMM Algorithm with the Uniform kernel 
+ * and a separation criterion of -1 (i.e. no nearfield and only farfield interaction)
+ * and compares the results with a direct computation.
+ * The matrix kernel has to be smooth at the origin, e.g. 1/r^2.
+ */
+
+// Simply create particles and try the kernels
+int main(int argc, char* argv[])