diff --git a/Examples/LagrangeInterpolationFMM.cpp b/Examples/LagrangeInterpolationFMM.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..74a0fdf8e313b165e6831449748dc273fd541b7b
--- /dev/null
+++ b/Examples/LagrangeInterpolationFMM.cpp
@@ -0,0 +1,219 @@
+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
+// 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".
+// ===================================================================================
+
+// ==== CMAKE =====
+// @FUSE_FFT
+// ==== Git =====
+// @SCALFMM_PRIVATE
+// ================
+
+#include <iostream>
+
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+
+#include "ScalFmmConfig.h"
+
+#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 "Containers/FOctree.hpp"
+
+#ifdef _OPENMP
+#include "Core/FFmmAlgorithmThread.hpp"
+#else
+#include "Core/FFmmAlgorithm.hpp"
+#endif
+
+#include "Utils/FParameterNames.hpp"
+
+/**
+ * This program runs the FMM Algorithm with the uniform interpolation kernel and compares the results with a direct computation.
+ */
+/// \file LagrangeInterpolationFMM.cpp
+//!
+//! \brief This program runs the FMM Algorithm with the interpolation kernel based on uniform (grid points) interpolation (1/r kernel)
+//! \authors B. Bramas, O. Coulaud
+//!
+//! This code is a short example to use the Interpolation approach for the 1/r kernel
+
+
+// Simply create particles and try the kernels
+int main(int argc, char* argv[])
+{
+ FHelpDescribeAndExit(argc, argv,
+ "Driver for Lagrange interpolation kernel (1/r kernel).",
+ FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
+ FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
+ FParameterDefinitions::NbThreads);
+
+
+ const std::string defaultFile(/*SCALFMMDataPath+*/"Data/unitCubeXYZQ100.bfma" );
+ const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
+ const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
+ const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
+ const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
+
+
+#ifdef _OPENMP
+ omp_set_num_threads(NbThreads);
+ std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
+#else
+ std::cout << "\n>> Sequential version.\n" << std::endl;
+#endif
+ //
+ std::cout << "Parameters "<< std::endl
+ << " Octree Depth "<< TreeHeight <<std::endl
+ << " SubOctree depth " << SubTreeHeight <<std::endl
+ << " Input file name: " <<filename <<std::endl
+ << " Thread number: " << NbThreads <<std::endl
+ <<std::endl;
+ //
+ // init timer
+ FTic time;
+
+ // open particle file
+ ////////////////////////////////////////////////////////////////////
+ //
+ FFmaGenericLoader loader(filename);
+ //
+ ////////////////////////////////////////////////////////////////////
+ // begin Lagrange kernel
+
+ // accuracy
+ const unsigned int ORDER = 7;
+ // typedefs
+ typedef FP2PParticleContainerIndexed<> ContainerClass;
+ typedef FSimpleLeaf< ContainerClass > LeafClass;
+ typedef FUnifCell<ORDER> CellClass;
+ typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
+ //
+ typedef FInterpMatrixKernelR MatrixKernelClass;
+ const MatrixKernelClass MatrixKernel;
+ typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+ //
+#ifdef _OPENMP
+ typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+#else
+ typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+#endif
+ // init oct-tree
+ OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+
+
+ { // -----------------------------------------------------
+ std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
+ << " particles ..." << std::endl;
+ std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
+ time.tic();
+ //
+ FPoint position;
+ FReal physicalValue = 0.0;
+ //
+ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+ //
+ // Read particle per particle from file
+ loader.fillParticle(&position,&physicalValue);
+ //
+ // put particle in octree
+ tree.insert(position, idxPart, physicalValue);
+ }
+
+ time.tac();
+ std::cout << "Done " << "(@Creating and Inserting Particles = "
+ << time.elapsed() << " s) ." << std::endl;
+ } // -----------------------------------------------------
+
+ { // -----------------------------------------------------
+ std::cout << "\nLagrange FMM (ORDER="<< ORDER << ") ... " << std::endl;
+
+ time.tic();
+ //
+ KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
+ //
+ FmmClass algo(&tree, &kernels);
+ //
+ algo.execute(); // Here the call of the FMM algorithm
+ //
+ time.tac();
+ std::cout << "Timers Far Field \n"
+ << "P2M " << algo.getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
+ << "M2M " << algo.getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
+ << "M2L " << algo.getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
+ << "L2L " << algo.getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
+ << "P2P and L2P " << algo.getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
+ << std::endl;
+
+
+ std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
+ }
+ // -----------------------------------------------------
+ //
+ // Some output
+ //
+ //
+ { // -----------------------------------------------------
+ long int N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
+ FReal energy =0.0 ;
+ //
+ // Loop over all leaves
+ //
+ std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
+ std::cout << std::scientific;
+ std::cout.precision(10) ;
+
+ tree.forEachLeaf([&](LeafClass* leaf){
+ const FReal*const posX = leaf->getTargets()->getPositions()[0];
+ const FReal*const posY = leaf->getTargets()->getPositions()[1];
+ const FReal*const posZ = leaf->getTargets()->getPositions()[2];
+
+ const FReal*const potentials = leaf->getTargets()->getPotentials();
+ const FReal*const forcesX = leaf->getTargets()->getForcesX();
+ const FReal*const forcesY = leaf->getTargets()->getForcesY();
+ const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
+ const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
+ const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
+
+ const FVector<int>& indexes = leaf->getTargets()->getIndexes();
+
+ for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+ const int indexPartOrig = indexes[idxPart];
+ if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
+ std::cout << "Index "<< indexPartOrig <<" potential " << potentials[idxPart]
+ << " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
+ << " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
+ }
+ energy += potentials[idxPart]*physicalValues[idxPart] ;
+ }
+ });
+ std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
+ std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
+
+ }
+ // -----------------------------------------------------
+
+
+ return 0;
+}