Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

Addons/CKernelApi/Src/FScalfmmApiInit.cpp

@@ -6,48 +6,94 @@ extern "C" {
 #include "FInterEngine.hpp"
 #include "FUserKernelEngine.hpp"
 
-extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double* BoxCenter, */scalfmm_kernel_type KernelType){
-    ScalFmmCoreHandle * handle = new ScalFmmCoreHandle();
+extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double* BoxCenter, */scalfmm_kernel_type KernelType,
+                                       scalfmm_algorithm algo){
+    ScalFmmCoreHandle<double> * handle = new ScalFmmCoreHandle<double>();
     typedef double FReal;
 
-    switch(KernelType){
-    case 0:
-        handle->engine = new FUserKernelEngine<FReal>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
-        break;
+    if(algo == source_target){
 
-    case 1:
-        //TODO typedefs
-        typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
-        typedef FChebCell<FReal,7>                                                         ChebCell;
+        switch(KernelType){
+        case 0:
+            typedef FP2PParticleContainerIndexed<FReal>           ContainerClass;
+            typedef FTypedLeaf<FReal,ContainerClass>                                         LeafClass;
 
-        typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
-        typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
+            handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
+            break;
 
-        handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
-        break;
-        // case 2:
-        //     //TODO typedefs
-        //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
-        //     typedef FUnifCell<7>                                                         UnifCell;
+        case 1:
+            //TODO typedefs
+            typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
+            typedef FTypedChebCell<FReal,7>                                                   ChebCell;
+            typedef FTypedLeaf<FReal,ContainerClass>                                         LeafClass;
 
-        //     typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
-        //     typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7>           UnifKernel;
+            typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
+            typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
 
-        //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
-        //     break;
+            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            break;
+            // case 2:
+            //     //TODO typedefs
+            //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
+            //     typedef FUnifCell<7>                                                         UnifCell;
 
-    default:
-        std::cout<< "Kernel type unsupported" << std::endl;
-        exit(0);
-        break;
+            //     typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
+            //     typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7>           UnifKernel;
+
+            //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            //     break;
+
+        default:
+            std::cout<< "Kernel type unsupported" << std::endl;
+            exit(0);
+            break;
+        }
+    }
+    else{ //No Source/Targets distinction
+        switch(KernelType){
+        case 0:
+            typedef FP2PParticleContainerIndexed<FReal>                            ContainerClass;
+            typedef FSimpleLeaf<FReal,ContainerClass>                                   LeafClass;
+
+            handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
+            break;
+
+        case 1:
+            //TODO typedefs
+            typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
+            typedef FChebCell<FReal,7>                                                   ChebCell;
+            typedef FSimpleLeaf<FReal,ContainerClass>                                         LeafClass;
+
+            typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
+            typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
+
+            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            break;
+            // case 2:
+            //     //TODO typedefs
+            //     typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
+            //     typedef FUnifCell<7>                                                         UnifCell;
+
+            //     typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
+            //     typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7>           UnifKernel;
+
+            //     handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            //     break;
+
+        default:
+            std::cout<< "Kernel type unsupported" << std::endl;
+            exit(0);
+            break;
+
+        }
     }
-    return handle;
+   return handle;
 }
 
 extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
-    ((ScalFmmCoreHandle *) handle)->engine->intern_dealloc_handle(userDeallocator);
-    delete ((ScalFmmCoreHandle *) handle)->engine ;
-    delete (ScalFmmCoreHandle *) handle;
+    ((ScalFmmCoreHandle<double> *) handle)->engine->intern_dealloc_handle(userDeallocator);
+    delete ((ScalFmmCoreHandle<double> *) handle)->engine ;
+    delete (ScalFmmCoreHandle<double> *) handle;
 }
 
 /**

Addons/CKernelApi/Tests/testSphereElectro.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+
+
+
+//For timing monitoring
+#include <time.h>
+#include <sys/time.h>
+
+#include "../Src/CScalfmmApi.h"
+
+#include "../../Src/Kernels/Chebyshev/FChebInterface.h"
+
+double getRandom(){
+    return (random()/(double)(RAND_MAX));
+}
+
+void generateSurfacePointOnUnitSphere(int N, double * points){
+    double u, v, theta, phi, sinPhi ;
+    //
+    int j = 0,i=0 ;
+    for ( i = 0 ; i< N ; ++i, j+=3)  {
+        //
+        u = getRandom() ;  v = getRandom() ;
+        theta  = 2*M_PI*u ;
+        phi     = acos(2*v-1);
+        sinPhi = sin(phi);
+        //
+        points[j]   =     cos(theta)*sinPhi ;
+        points[j+1] =   sin(theta)*sinPhi ;
+        points[j+2] =   2*v-1 ;
+        //
+    }
+}
+
+void generateSurfacePoints(double rayon, double* centre, int nbDePoints, double* points){
+
+    generateSurfacePointOnUnitSphere(nbDePoints , points) ;
+    int j =0,i=0 ;
+    for ( i = 0 ; i< nbDePoints ; ++i, j+=3)  {
+        points[j]    *= rayon + centre[0];
+        points[j+1]  *= rayon + centre[1];
+        points[j+2]  *= rayon + centre[2];
+    }
+}
+
+void generateInsidePoints(double rayon, double*centre, int nbDePoints,  double* points){
+    generateSurfacePointOnUnitSphere(nbDePoints, points);
+    int j=0;
+    double u;
+    for(j=0 ; j<nbDePoints ; ++j){
+        u = getRandom();
+        points[j]    *= (rayon + centre[0])*u;
+        points[j+1]  *= (rayon + centre[1])*u;
+        points[j+2]  *= (rayon + centre[2])*u;
+    }
+}
+
+void displayPoints(int nbPoints, double * points){
+    int i = 0;
+    for(i=0 ; i<nbPoints ; ++i){
+        printf("%e %e %e \n",points[i*3],points[i*3+1],points[i*3+2]);
+    }
+}
+
+void displayArray(int nbValue, double * array){
+    int i = 0;
+    for(i=0 ; i<nbValue ; ++i){
+        printf("%e \n",array[i]);
+    }
+}
+
+void getNormal(double * positions, double * normeToFill){
+    int i;
+    double norme = sqrt(positions[0]*positions[0] + positions[1]*positions[1] + positions[2]*positions[2]);
+    for(i=0 ; i<3 ; ++i){
+        normeToFill[i] = positions[i]/norme;
+    }
+    printf("Tgt Norme %e - %e - %e\n",
+           normeToFill[0],
+           normeToFill[1],
+           normeToFill[2]);
+}
+
+void computeNormalXForces(int nbPoints, double * forcesToRead, double * positionsToRead, double * arrayToFill){
+    double * currentNormal = malloc(sizeof(double)*3);
+    int idxPart,i;
+    for(idxPart = 0 ; idxPart<nbPoints ; ++idxPart){
+        getNormal(&positionsToRead[idxPart],currentNormal); //get the norme
+        for(i=0 ; i<3 ; ++i){
+            arrayToFill[idxPart] += currentNormal[i]*forcesToRead[idxPart+i];
+        }
+    }
+    free(currentNormal);
+}
+
+int main(int argc, char ** av){
+    printf("Start\n");
+    if(argc<2){
+        printf("Use : %s nb_part(cible) (optionnal : TreeHeight) \nexiting\n",av[0]);
+        exit(0);
+    }
+    int nbPartTarget= atoi(av[1]);
+    int treeHeight = 5 ;
+    if(argc>2){
+        int treeHeight = atoi(av[2]);
+    }
+
+    double boxWidth = 2.0;
+    double boxCenter[3];
+    boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0;
+
+    int i;
+    //Allocation of the target points
+    double * targetsXYZ = malloc(sizeof(double)* 3*nbPartTarget);
+    double * targetsPhiValues = malloc(sizeof(double)* nbPartTarget);
+    //Memset (au cas ou)
+    memset(targetsXYZ,0,sizeof(double)*3*nbPartTarget);
+    memset(targetsPhiValues,0,sizeof(double)*nbPartTarget);
+    //Fill
+    for(i=0 ; i<nbPartTarget ; ++i){
+        targetsPhiValues[i] = -1.0;
+    }
+    generateSurfacePoints(1.0,boxCenter,nbPartTarget,targetsXYZ);
+    printf("Surface points generated \n");
+
+    //Allocation of the sources points
+    int nbPartSource = 10;
+    double * sourceXYZ = malloc(sizeof(double)* 3*nbPartSource);
+    double * sourcePhiValues = malloc(sizeof(double)* nbPartSource);
+    //Set to Zero
+    memset(sourceXYZ,0,3*sizeof(double)*nbPartSource);
+    memset(sourcePhiValues,0,sizeof(double)*nbPartSource);
+    //Fill
+    for(i=0 ; i<nbPartSource ; ++i){
+        sourcePhiValues[i] = 1.0;
+    }
+    generateInsidePoints(1.0,boxCenter,nbPartSource,sourceXYZ);
+    //displayPoints(nbPartTarget,targetsXYZ);
+
+    printf("Inside points generated \n");
+    //displayPoints(nbPartSource,sourceXYZ);
+    //Creation of arrays to store forces
+    double * arrayOfForces = malloc(sizeof(double )* 3 * (nbPartSource+nbPartTarget));
+    memset(arrayOfForces,0,sizeof(double)* 3 * (nbPartTarget));
+
+    {//Start of computation
+
+        //For handling the library
+        scalfmm_handle handle = scalfmm_init(chebyshev,source_target);
+
+        //Struct for ref cheb kernel
+        struct User_Scalfmm_Cell_Descriptor user_descr;
+        user_descr.user_init_cell = NULL;
+        user_descr.user_free_cell = NULL;
+        //Set algorithm to source target
+        //scalfmm_algorithm_config(handle,source_target);
+        //Build the tree
+        scalfmm_build_tree(handle,treeHeight, boxWidth, boxCenter, user_descr);
+
+        //Insert Sources and targets
+        scalfmm_tree_insert_particles_xyz(handle,nbPartSource,sourceXYZ,SOURCE);
+        printf("Sources inserted \n");
+        scalfmm_tree_insert_particles_xyz(handle,nbPartTarget,targetsXYZ,TARGET);
+        printf("Targets inserted \n");
+        //Since we inserted first the sources, then sources will get
+        //indices from 0 to (nbPartSource-1), and targets from
+        //(nbPartSource) to nbPartSource+nbPartTarget-1).
+
+        int * arrayofIndicesSource = malloc(sizeof(int)*nbPartSource);
+        int * arrayofIndicesTarget = malloc(sizeof(int)*nbPartTarget);
+        {//Set physical values
+
+            //SRC
+            int idPart;
+            for(idPart = 0 ; idPart<nbPartSource ; ++idPart){
+                arrayofIndicesSource[idPart] = idPart;
+            }
+            scalfmm_set_physical_values_npart(handle,nbPartSource,arrayofIndicesSource,sourcePhiValues,SOURCE);
+            //TGT
+            for(idPart = 0 ; idPart<nbPartTarget ; ++idPart){
+                arrayofIndicesTarget[idPart] = idPart; // here, we add the number of sources previously inserted
+            }
+            scalfmm_set_physical_values_npart(handle,nbPartTarget,arrayofIndicesTarget,targetsPhiValues,TARGET);
+
+
+
+        }
+        //Computation
+        scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
+
+        //Get back the forces
+        scalfmm_get_forces_xyz(handle,nbPartTarget,arrayOfForces,TARGET);
+        scalfmm_get_forces_xyz(handle,nbPartSource,&arrayOfForces[nbPartTarget],SOURCE);
+        printf("Forces computed : \n");
+        displayPoints(nbPartTarget+nbPartSource,arrayOfForces);
+        printf("As expected, Source forces are 0\n \n");
+        //Release memory used :
+        free(arrayofIndicesSource);
+        free(arrayofIndicesTarget);
+
+        scalfmm_dealloc_handle(handle,NULL);
+
+    }
+
+    {//Let's check the result, we computed fr each target part its forces
+        //Storage of reference forces
+        double * arrayRefForces = malloc(sizeof(double)*nbPartTarget*3);
+        memset(arrayRefForces,0,sizeof(double)*nbPartTarget*3);
+
+        int idTgt;
+        for(idTgt = 0 ; idTgt<nbPartTarget ; ++idTgt){
+            int idSrc;
+            double dx,dy,dz;
+            for(idSrc = 0 ; idSrc<nbPartTarget ; ++idSrc){
+                //First compute dist.
+                dx = sourceXYZ[idSrc+0] - targetsXYZ[idTgt+0];
+                dy = sourceXYZ[idSrc+1] - targetsXYZ[idTgt+1];
+                dz = sourceXYZ[idSrc+2] - targetsXYZ[idTgt+2];
+
+                //Secondly, compute coeff
+                double coeffs = targetsPhiValues[idTgt] * sourcePhiValues[idSrc];
+                double one_over_r = 1.0/(sqrt(dx*dx+dy*dy+dz*dz));
+                double one_over_r3 = one_over_r * one_over_r * one_over_r;
+
+                arrayRefForces[idTgt*3+0] += dx*coeffs*one_over_r3;
+                arrayRefForces[idTgt*3+1] += dy*coeffs*one_over_r3;
+                arrayRefForces[idTgt*3+2] += dz*coeffs*one_over_r3;
+
+            }
+        }
+
+        {//Then, we compare
+            double errorCumul = 0;
+            int idArr;
+            for(idArr = 0 ; idArr<nbPartTarget ; ++idArr){
+                errorCumul += fabs(arrayRefForces[idArr+0]-arrayOfForces[idArr+0]);
+                errorCumul += fabs(arrayRefForces[idArr+1]-arrayOfForces[idArr+1]);
+                errorCumul += fabs(arrayRefForces[idArr+2]-arrayOfForces[idArr+2]);
+                printf("Directly Computed %e %e %e\n",
+                       arrayRefForces[idArr+0],
+                       arrayRefForces[idArr+1],
+                       arrayRefForces[idArr+2]);
+            }
+            printf("Error cumul : %e\n",errorCumul);
+        }
+    }
+
+
+    //Part where we apply normal on target's forces vector
+    //Copying each target's parts forces,
+    double * targetsForces = malloc(sizeof(double) * 3 * nbPartTarget);
+    memcpy(targetsForces,arrayOfForces,sizeof(double)*3*nbPartTarget);
+
+    double * normeXForces =  malloc(sizeof(double) * nbPartTarget);
+    memset(normeXForces,0,sizeof(double) * nbPartTarget);
+
+    computeNormalXForces(nbPartTarget,targetsForces,targetsXYZ,normeXForces);
+    printf("For each target, we display [Normal To Sphere] . [Force product] \n");
+    displayArray(nbPartTarget,normeXForces);
+
+
+    //Free memory
+    free(sourceXYZ);
+    free(sourcePhiValues);
+    free(targetsXYZ);
+    free(targetsPhiValues);
+    free(arrayOfForces);
+    return EXIT_SUCCESS;
+}

CMakeModules/morse/find/FindBLAS.cmake

@@ -301,6 +301,10 @@ macro(Check_Fortran_Libraries LIBRARIES _prefix _name _flags _list _thread)
             unset(${_prefix}${_combined_name}_WORKS CACHE)
         endif()
         if (_CHECK_FORTRAN)
+            if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+                string(REPLACE "mkl_intel_lp64" "mkl_gf_lp64" CMAKE_REQUIRED_LIBRARIES "${CMAKE_REQUIRED_LIBRARIES}")
+                string(REPLACE "mkl_intel_ilp64" "mkl_gf_ilp64" CMAKE_REQUIRED_LIBRARIES "${CMAKE_REQUIRED_LIBRARIES}")
+            endif()
             check_fortran_function_exists("${_name}" ${_prefix}${_combined_name}_WORKS)
         else()
             check_function_exists("${_name}_" ${_prefix}${_combined_name}_WORKS)
@@ -415,8 +419,6 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
         endif()
     endif()
 
-    # libiomp5
-    # --------
     if (WIN32)
         string(REPLACE ":" ";" _libdir "$ENV{LIB}")
     elseif (APPLE)
@@ -426,6 +428,8 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
     endif ()
     list(APPEND _libdir "${CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES}")
     list(APPEND _libdir "${CMAKE_C_IMPLICIT_LINK_DIRECTORIES}")
+    # libiomp5
+    # --------
     set(OMP_iomp5_LIBRARY "OMP_iomp5_LIBRARY-NOTFOUND")
     find_library(OMP_iomp5_LIBRARY
         NAMES iomp5
@@ -433,15 +437,35 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
       )
     mark_as_advanced(OMP_iomp5_LIBRARY)
     set(OMP_LIB "")
-    if (OMP_iomp5_LIBRARY)
-        set(OMP_LIB "${OMP_iomp5_LIBRARY}")
+    # libgomp
+    # -------
+    set(OMP_gomp_LIBRARY "OMP_gomp_LIBRARY-NOTFOUND")
+    find_library(OMP_gomp_LIBRARY
+        NAMES gomp
+        HINTS ${_libdir}
+      )
+    mark_as_advanced(OMP_gomp_LIBRARY)
+    # choose one or another depending on the compilo
+    if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
+        if (OMP_gomp_LIBRARY)
+            set(OMP_LIB "${OMP_gomp_LIBRARY}")
+        endif()
+    else(CMAKE_C_COMPILER_ID STREQUAL "Intel")
+        if (OMP_iomp5_LIBRARY)
+            set(OMP_LIB "${OMP_iomp5_LIBRARY}")
+        endif()
     endif()
 
     if (UNIX AND NOT WIN32)
         set(LM "-lm")
         set(BLAS_COMPILER_FLAGS "")
-        if (CMAKE_C_COMPILER_ID STREQUAL "Intel" AND NOT BLA_VENDOR STREQUAL "Intel10_64lp_seq")
-            list(APPEND BLAS_COMPILER_FLAGS "-openmp")
+        if (NOT BLA_VENDOR STREQUAL "Intel10_64lp_seq")
+            if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
+                list(APPEND BLAS_COMPILER_FLAGS "-openmp")
+            endif()
+            if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
+                list(APPEND BLAS_COMPILER_FLAGS "-fopenmp")
+            endif()
         endif()
         if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
             if (BLA_VENDOR STREQUAL "Intel10_32")
@@ -523,8 +547,14 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
                     list(APPEND BLAS_SEARCH_LIBS
                     "mkl_blas95 mkl_intel_lp64 mkl_intel_thread mkl_core guide")
                     # mkl >= 10.3
-                    list(APPEND BLAS_SEARCH_LIBS
-                    "mkl_blas95_lp64 mkl_intel_lp64 mkl_intel_thread mkl_core")
+                    if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
+                        list(APPEND BLAS_SEARCH_LIBS
+                        "mkl_blas95_lp64 mkl_intel_lp64 mkl_intel_thread mkl_core")
+                    endif()
+                    if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
+                        list(APPEND BLAS_SEARCH_LIBS
+                        "mkl_blas95_lp64 mkl_intel_lp64 mkl_gnu_thread mkl_core")
+                    endif()
                 endif ()
                 if (BLA_VENDOR STREQUAL "Intel10_64lp_seq" OR BLA_VENDOR STREQUAL "All")
                     list(APPEND BLAS_SEARCH_LIBS
@@ -589,8 +619,14 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
                     list(APPEND BLAS_SEARCH_LIBS
                     "mkl_intel_lp64 mkl_intel_thread mkl_core guide")
                     # mkl >= 10.3
-                    list(APPEND BLAS_SEARCH_LIBS
-                    "mkl_intel_lp64 mkl_intel_thread mkl_core")
+                    if (CMAKE_C_COMPILER_ID STREQUAL "Intel")
+                        list(APPEND BLAS_SEARCH_LIBS
+                        "mkl_intel_lp64 mkl_intel_thread mkl_core")
+                    endif()
+                    if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
+                        list(APPEND BLAS_SEARCH_LIBS
+                        "mkl_intel_lp64 mkl_gnu_thread mkl_core")
+                    endif()
                 endif ()
                 if (BLA_VENDOR STREQUAL "Intel10_64lp_seq" OR BLA_VENDOR STREQUAL "All")
                     list(APPEND BLAS_SEARCH_LIBS

CMakeModules/morse/find/FindMETIS.cmake

@@ -205,6 +205,10 @@ if(METIS_LIBRARIES)
         set(REQUIRED_LIBDIRS "${METIS_LIBRARY_DIRS}")
     endif()
     set(REQUIRED_LIBS "${METIS_LIBRARIES}")
+    # m
+    if(UNIX OR WIN32)
+        list(APPEND REQUIRED_LIBS "-lm")
+    endif()
 
     # set required libraries for link
     set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")

CMakeModules/morse/find/FindPARMETIS.cmake

@@ -205,6 +205,10 @@ if(PARMETIS_LIBRARIES)
         set(REQUIRED_LIBDIRS "${PARMETIS_LIBRARY_DIRS}")
     endif()
     set(REQUIRED_LIBS "${PARMETIS_LIBRARIES}")
+    # m
+    if(UNIX OR WIN32)
+        list(APPEND REQUIRED_LIBS "-lm")
+    endif()
 
     # set required libraries for link
     set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")

Src/Components/FTypedLeaf.hpp

@@ -4,13 +4,13 @@
 // 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.  
-// 
+// 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.cecill.info".
 // "http://www.gnu.org/licenses".
 // ===================================================================================
 #ifndef FTYPEDLEAF_HPP
@@ -55,6 +55,20 @@ public:
         else sources.push(inParticlePosition, FParticleTypeSource, args...);
     }
 
+    /**
+     * To add a new particle in the leaf
+     * @param inParticlePosition the position of the new particle
+     * @param isTarget bool to know if it is a target
+     * followed by other param given by the user
+     */
+    template<typename... Args>
+    void push(const FPoint<FReal>& inParticlePosition, Args ... args){
+        FAssert(0,"Error : cannot push a particle without specifying type (src/tgt)");
+    }
+
+
+
+
     /**
     * To get all the sources in a leaf
     * @return a pointer to the list of particles that are sources
@@ -75,5 +89,3 @@ public:
 
 
 #endif //FTYPEDLEAF_HPP
-
-

Src/Containers/FBoolArray.hpp

@@ -16,7 +16,8 @@
 #ifndef FBOOLARRAY_HPP
 #define FBOOLARRAY_HPP
 
-
+#include "../Utils/FGlobal.hpp"
+#include "../Utils/FAssert.hpp"
 // To get memcpy
 #include <cstring>
 
@@ -31,29 +32,29 @@
 */
 class FBoolArray{
     /** Size of a unsigned long */
-    const static int BytesInBlock = sizeof(unsigned long);
-    const static int SizeOfBlock = BytesInBlock * 8;
+    const static FSize BytesInBlock = sizeof(unsigned long);
+    const static FSize SizeOfBlock = BytesInBlock * 8;
 
     /** The array to store bits */
-    unsigned long* const array;
+    unsigned long* array;
     /** Size of the memory allocated */
-    const int memSize;
+    FSize memSize;
     /** Size of the array => number of real elements */
-    const int size;
+    FSize size;
 
     /** get size to number of long */
-    int LongFromSize(const int inSize){
+    FSize LongFromSize(const FSize inSize){
         return ((inSize + SizeOfBlock - 1) / SizeOfBlock);
     }
 
     /** Alloc an array */
-    unsigned long * AllocArray(const int inSize){
+    unsigned long * AllocArray(const FSize inSize){
         return new unsigned long[LongFromSize(inSize)];
     }
 
 public :
     /** Constructor with size */
-    FBoolArray(const int inSize) : array(AllocArray(inSize)), memSize(LongFromSize(inSize)*BytesInBlock), size(inSize) {
+    explicit FBoolArray(const FSize inSize = 0) : array(AllocArray(inSize)), memSize(LongFromSize(inSize)*BytesInBlock), size(inSize) {
         setToZeros();
     }
 
@@ -62,6 +63,25 @@ public :
         *this = other;
     }
 
+    /** Move the data */
+    FBoolArray(FBoolArray&& other): array(nullptr), memSize(0), size(0){
+        array   = other.array;
+        memSize = other.memSize;
+        size    = other.size;
+        other.array   = nullptr;
+        other.memSize = 0;
+        other.size    = 0;
+    }
+
+    /** remove all values and allocate new array */
+    void reset(const FSize inSize){
+        delete [] array;
+        array   = (AllocArray(inSize));
+        memSize = (LongFromSize(inSize)*BytesInBlock);
+        size    = (inSize);
+        setToZeros();
+    }
+
     /** Destructor */
     ~FBoolArray(){
         delete [] array;
@@ -72,10 +92,25 @@ public :
     * Array must have the same size
     */
     FBoolArray& operator=(const FBoolArray& other){
+        FAssertLF(size == other.size);
         memcpy(array, other.array, memSize);
         return *this;
     }
 
+    /**
+     * Move the data from one array to the other
+     */
+    FBoolArray& operator=(FBoolArray&& other){
+        delete [] array;
+        array   = other.array;
+        memSize = other.memSize;
+        size    = other.size;
+        other.array   = nullptr;
+        other.memSize = 0;
+        other.size    = 0;
+        return *this;
+    }
+
     /**
     * Operator ==
     * Array must have the same size
@@ -93,22 +128,22 @@ public :
     }
 
     /** To get a value */
-    bool get(const int inPos) const {
-        const int posInArray = inPos / SizeOfBlock;
-        const int bytePosition = inPos - (posInArray * 8);
+    bool get(const FSize inPos) const {
+        const FSize posInArray = inPos / SizeOfBlock;
+        const FSize bytePosition = inPos - (posInArray * 8);
         return (array[posInArray] >> bytePosition) & 1;
     }
 
     /** To set a value */
-    void set(const int inPos, const bool inVal){
-        const int posInArray = inPos / SizeOfBlock;
-        const int bytePosition = inPos - (posInArray * 8);
+    void set(const FSize inPos, const bool inVal){
+        const FSize posInArray = inPos / SizeOfBlock;
+        const FSize bytePosition = inPos - (posInArray * 8);
         if(inVal) array[posInArray] |= (1UL << bytePosition);
         else array[posInArray] &= ~(1UL << bytePosition);
     }
 
     /** To get the size of the array */
-    int getSize() const {