New feature : It's possible to end the FMM at a given level greater than 2. If...

New feature : It's possible to end the FMM at a given level greater than 2. If used, a specific M2L operator (given by user) will be applied on that last given level between every pairs of cells

Addons/CKernelApi/Src/CScalfmmApi.h

@@ -429,6 +429,17 @@ typedef void (*Callback_M2M)(int level, void* parentCell, int childPosition, voi
  */
 typedef void (*Callback_M2L)(int level, void* targetCell, int sourceCellPosition, void* sourceCell, void* userData);
 
+/**
+ * @brief Function to be filled by user's M2L
+ * @param level current level in the tree
+ * @param targetCell pointer to cell to be filled
+ * @param sourceCell cell to be read
+ * @param transfer array of 3 int, displaying the number of boxes at
+ * the current level in along X,Y and Z axis.
+ * @param userData datas specific to the user's kernel
+ */
+typedef void (*Callback_M2L_Ext)(int level, void* targetCell, void* sourceCell, int tranfer[3], void* userData);
+
 /**
  * @brief Function to be filled by user's M2L
  * @param level current level in the tree
@@ -508,6 +519,7 @@ typedef struct User_Scalfmm_Kernel_Descriptor {
     Callback_P2M p2m;
     Callback_M2M m2m;
     Callback_M2L m2l;
+    Callback_M2L_Ext m2l_ext;
     Callback_M2LFull m2l_full;
     Callback_L2L l2l;
     Callback_L2P l2p;
@@ -600,5 +612,14 @@ int scalfmm_get_nb_timers(scalfmm_handle handle);
 void scalfmm_get_timers(scalfmm_handle handle,double * Timers);
 
 
+/**
+ * @brief Set the upper limit int the tree for applying FMM : standard
+ * value is 2. If used, then user MUST provide a M2L source to target
+ * to be applyied on a wider range than the usual 189 cells.
+ * @param Handle scalfmm_handle provided by scalfmm_init.
+ * @param upperLimit : int  than 2.
+ */
+void scalfmm_set_upper_limit(scalfmm_handle handle, int upperLimit);
+
 
 #endif

Addons/CKernelApi/Src/FScalFMMEngine.hpp

@@ -824,6 +824,11 @@ public:
     virtual int get_nb_timers(){
         return algoTimer->getNbOfTimerRecorded();
     }
+
+    virtual void set_upper_limit(int upperLimit){
+        FAssertLF(0,"This feature is not available with Chebyshev Kernel, please use your own kernel or do not use it.\n Exiting anyways...\n");
+    }
+
 };
 
 template<class FReal>
@@ -1045,4 +1050,8 @@ extern "C" void scalfmm_print_everything(scalfmm_handle Handle){
     ((ScalFmmCoreHandle<double> * ) Handle)->engine->print_everything();
 }
 
+extern "C" void scalfmm_set_upper_limit(scalfmm_handle Handle, int upperLimit){
+    ((ScalFmmCoreHandle<double> * ) Handle)->engine->set_upper_limit(upperLimit);
+}
+
 #endif

Addons/CKernelApi/Src/FScalfmmApiInit.cpp

@@ -103,6 +103,10 @@ extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell
  */
 #ifndef CHEBINTERFACE_HPP
 #define CHEBINTERFACE_HPP
+
+#warning "Compiling Cheb Interface"
+
+
 extern "C" {
 #include "Kernels/Chebyshev/FChebInterface.h"
 }

Addons/CKernelApi/Src/FUserKernelEngine.hpp

@@ -96,6 +96,7 @@ public:
 
     /** Default destructor */
     virtual ~CoreKernel(){
+
     }
 
     /** Do nothing */
@@ -108,7 +109,6 @@ public:
         if(kernel.m2m){
             for(int idx = 0 ; idx < 8 ; ++idx){
                 if( children[idx] ){
-                    printf("lvl : %d\n",level);
                     kernel.m2m(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
                 }
             }
@@ -201,6 +201,13 @@ public:
         return userData;
     }
 
+    void M2L_Extended(CellClass * src, CellClass * tgt, const FTreeCoordinate transfer, const int level){
+        if(kernel.m2l_ext){
+            int array[3] = {transfer.getX(),transfer.getY(),transfer.getZ()};
+            kernel.m2l_ext(level,tgt->getContainer(),src->getContainer(),array,userData);
+        }
+    }
+
 };
 
 template<class FReal,class LeafClass>
@@ -212,7 +219,6 @@ private:
 
     //Typedefs :
     typedef FOctree<FReal,CoreCell,ContainerClass,LeafClass>            OctreeClass;
-
     typedef CoreKernel<CoreCell,ContainerClass>           CoreKernelClass;
 
     //For arranger classes
@@ -220,6 +226,8 @@ private:
     //Attributes
     OctreeClass * octree;
     CoreKernelClass * kernel;
+    int upperLimit;
+    int treeHeight;
 
     // ArrangerClass * arranger;
     // ArrangerClassTyped * arrangerTyped;
@@ -228,7 +236,7 @@ private:
 
 public:
     FUserKernelEngine(/*int TreeHeight, double BoxWidth , double * BoxCenter, */scalfmm_kernel_type KernelType) :
-        octree(nullptr), kernel(nullptr) /*,arranger(nullptr)*/ {
+        octree(nullptr), kernel(nullptr), upperLimit(2), treeHeight(0) /*,arranger(nullptr)*/ {
         FScalFMMEngine<FReal>::kernelType = KernelType;
     }
 
@@ -254,6 +262,7 @@ public:
 
     void build_tree(int TreeHeight,double BoxWidth,double* BoxCenter,Scalfmm_Cell_Descriptor user_cell_descriptor){
         CoreCell::Init(user_cell_descriptor);
+        this->treeHeight = TreeHeight;
         printf("Tree Height : %d \n",TreeHeight);
         this->octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint<FReal>(BoxCenter));
     }
@@ -437,23 +446,106 @@ public:
             });
     }
 
+    void set_upper_limit(int inUpperLimit){
+        upperLimit = inUpperLimit;
+    }
+
+    /**
+     * @brief This function is called if the FMM is not computed on
+     * all the standards levels
+     *
+     */
+    void internal_M2L(){
+        if(upperLimit > 1){ // if upperLimit == 1, then, M2L has been
+                            // done at level 2, and hence all the far
+                            // field has been calculated.
+            //Starting at the lower level where the M2L has not been done.
+            typename OctreeClass::Iterator octreeIterator(octree); //lvl : 1
+
+            while(octreeIterator.level() != upperLimit){
+                octreeIterator.moveDown();
+            }
+
+            //I'm at the upperLimit, so the lowest level where M2L has been done.
+            do{
+                CoreCell * currentTgt = octreeIterator.getCurrentCell(); // This one is targeted
+
+                //Then, we get the interaction list at this lvl. This will provide us with lots of source cells.
+                const CoreCell * currentInteractionList[343];
+                //Get an iterator for the sources
+                typename OctreeClass::Iterator upAndDownIterator = octreeIterator;
+
+                {//This is supposed to be done for multiple level. You
+                 //need to go up until level 2. And then, to go down
+                 //until level upperLimit. I think it's possible ...
+                    while(upAndDownIterator.level() >= 2){
+                        upAndDownIterator.moveUp();
+
+                        //There, we get the interaction list of all parents of tgt cell
+                        const int nbInteract = octree->getInteractionNeighbors(currentInteractionList,
+                                                                               upAndDownIterator.getCurrentGlobalCoordinate(),
+                                                                               upAndDownIterator.level());
+                        int currentLevel = upAndDownIterator.level();
+                        if(nbInteract){
+                            //Then, we do M2L for each child at level upperLimit of each 343 Interaction cells.
+                            for(int idxSrc = 0; idxSrc < 343 ; ++idxSrc){
+                                if(currentInteractionList[idxSrc]){//Check if it exist
+                                    const CoreCell * currentSource = currentInteractionList[idxSrc]; //For clarity, will be otpimised out, anyway
+                                    MortonIndex idx = currentSource->getMortonIndex();
+
+                                    //At this point, we instanciate
+                                    //the number of child needed.
+                                    //This only depends on diffenrence
+                                    //between current level and
+                                    //upperLimit level
+                                    int totalNumberOfChild = FMath::pow(8,upperLimit-currentLevel);
+
+                                    for(int idxChildSrc = 0; idxChildSrc < totalNumberOfChild ; ++idxChildSrc){//For all 8^{number of levels to down} children
+                                        MortonIndex indexOfChild = ((idx << 3*(upperLimit-currentLevel))+idxChildSrc);
+                                        CoreCell * src = octree->getCell(indexOfChild,upperLimit); //Get the cell
+                                        if(src){//check if it exists
+                                            FTreeCoordinate srcCoord = src->getCoordinate();
+                                            FTreeCoordinate tgtCoord = currentTgt->getCoordinate();
+                                            //Build tree coord translation vector
+                                            FTreeCoordinate transfer;
+                                            transfer.setPosition(tgtCoord.getX()-srcCoord.getX(),
+                                                                 tgtCoord.getY()-srcCoord.getY(),
+                                                                 tgtCoord.getZ()-srcCoord.getZ());
+                                            kernel->M2L_Extended(src,currentTgt,transfer,octreeIterator.level());
+                                        }
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }while(octreeIterator.moveRight());
+        }
+        else{
+            FAssertLF("No reasons to be there, seriously ...\nExiting anyway...");
+        }
+    }
+
     void execute_fmm(){
         FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
+        FAbstractAlgorithm * abstrct = nullptr;
         switch(FScalFMMEngine<FReal>::Algorithm){
         case 0:
             {
                 typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassSeq;
                 AlgoClassSeq * algoSeq = new AlgoClassSeq(octree,kernel);
-                algoSeq->execute();
                 FScalFMMEngine<FReal>::algoTimer = algoSeq;
+                abstrct = algoSeq;
+                //algoSeq->execute(); will be done later
                 break;
             }
         case 1:
             {
                 typedef FFmmAlgorithmThread<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassThread;
                 AlgoClassThread*  algoThread = new AlgoClassThread(octree,kernel);
-                algoThread->execute();
                 FScalFMMEngine<FReal>::algoTimer = algoThread;
+                abstrct = algoThread;
+                //algoThread->execute(); will be done later
                 break;
             }
         case 2:
@@ -466,16 +558,30 @@ public:
             }
         case 3:
             {
-                // typedef FFmmAlgorithmThreadTsm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassTargetSource;
-                // AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
-                // algoTS->execute();
-                // FScalFMMEngine<FReal>::algoTimer = algoTS;
-                // break;
+                typedef FFmmAlgorithmThreadTsm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassTargetSource;
+                AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
+                FScalFMMEngine<FReal>::algoTimer = algoTS;
+                abstrct = algoTS;
+                //algoTS->execute(); will be done later
+                break;
             }
         default :
             std::cout<< "No algorithm found (probably for strange reasons) : "<< FScalFMMEngine<FReal>::Algorithm <<" exiting" << std::endl;
         }
-
+        if (FScalFMMEngine<FReal>::Algorithm != 2){
+            if(upperLimit != 2){
+                printf("At least I'm here\n");
+                abstrct->execute(FFmmP2M | FFmmM2M | FFmmM2L, upperLimit, treeHeight);
+                printf("\tUpPass finished\n");
+                internal_M2L();
+                printf("\tStrange M2L finished\n");
+                abstrct->execute(FFmmL2L | FFmmL2P | FFmmP2P, upperLimit, treeHeight);
+                printf("\tDownPass finished\n");
+            }
+            else{
+                abstrct->execute();
+            }
+        }
     }
 
     void intern_dealloc_handle(Callback_free_cell userDeallocator){

Addons/CKernelApi/Tests/testChebInterface.c

@@ -116,7 +116,7 @@ int main(int argc, char ** av){
     }
 
     scalfmm_handle handle = scalfmm_init(user_defined_kernel,multi_thread);
-
+    scalfmm_set_upper_limit(handle,2);
     //For Reference
     scalfmm_handle handle_ref = scalfmm_init(chebyshev,multi_thread);
 
@@ -213,6 +213,17 @@ int main(int argc, char ** av){
         scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
         tac(&interface_timer);
 
+        { //Temporary
+            int nbTimers = scalfmm_get_nb_timers(handle);
+            double * timersArray = malloc(sizeof(double)*nbTimers);
+            scalfmm_get_timers(handle,timersArray);
+            int i;
+            for(i=0; i<nbTimers ; ++i){
+                printf("ScalFMM Operands : %d : \t %e\n",i,timersArray[i]);
+            }
+            free(timersArray);
+        }
+
 
         //Reduction on forces & potential arrays
         {
@@ -236,6 +247,18 @@ int main(int argc, char ** av){
         scalfmm_execute_fmm(handle_ref/*, kernel, &my_data*/);
         tac(&ref_timer);
 
+        { //Temporary
+            int nbTimers = scalfmm_get_nb_timers(handle_ref);
+            double * timersArray = malloc(sizeof(double)*nbTimers);
+            scalfmm_get_timers(handle_ref,timersArray);
+            int i;
+            for(i=0; i<nbTimers ; ++i){
+                printf("ScalFMM Operands : %d : \t %e\n",i,timersArray[i]);
+            }
+            free(timersArray);
+        }
+
+
         printf("Intern Chebyshev done\n");
         print_elapsed(&ref_timer);
 

Addons/CKernelApi/Tests/testSphereElectro.c

@@ -178,7 +178,9 @@ int main(int argc, char ** av){
         int * arrayofIndicesSource = malloc(sizeof(int)*nbPartSource);
         int * arrayofIndicesTarget = malloc(sizeof(int)*nbPartTarget);
         {//Set physical values
-
+            printf("Setting physical values ... \n");
+            Timer phy_timer;
+            tic(&phy_timer);
             //SRC
             int idPart;
             for(idPart = 0 ; idPart<nbPartSource ; ++idPart){
@@ -190,15 +192,18 @@ int main(int argc, char ** av){
                 arrayofIndicesTarget[idPart] = idPart; // here, we add the number of sources previously inserted
             }
             scalfmm_set_physical_values_npart(handle,nbPartTarget,arrayofIndicesTarget,targetsPhiValues,TARGET);
-
-
+            tac(&phy_timer);
+            print_elapsed(&phy_timer);
+            printf("Physical values setted ... \n");
 
         }
+        printf("Executing FMM ...\n");
         Timer fmm_timer;
         tic(&fmm_timer);
         //Computation
         scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
         tac(&fmm_timer);
+        printf("FMM finished ...\n");
         print_elapsed(&fmm_timer);
 
         //Get back the forces