diff --git a/CMakeLists.txt b/CMakeLists.txt index 3ad1ec953ca421f8fb05efd101708b208978431d..0569a9306a16298d72b67f15efd44e4a486c8440 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ STRING(COMPARE EQUAL "${CMAKE_SOURCE_DIR}" "${CMAKE_BINARY_DIR}" insource) if(insource) MESSAGE(FATAL_ERROR "${PROJECT_NAME} requires an out of source build. Goto ./Build and tapes cmake ../") endif(insource) -SET(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/Modules) +SET(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMakeModules) # # MPI option has to be set before project, cannot be changed in the cache! diff --git a/Modules/CodeCoverage.cmake b/CMakeModules/CodeCoverage.cmake similarity index 100% rename from Modules/CodeCoverage.cmake rename to CMakeModules/CodeCoverage.cmake diff --git a/Modules/GetCpuInfos.cmake b/CMakeModules/GetCpuInfos.cmake similarity index 100% rename from Modules/GetCpuInfos.cmake rename to CMakeModules/GetCpuInfos.cmake diff --git a/Modules/checkAVXpe.cpp b/CMakeModules/checkAVXpe.cpp similarity index 100% rename from Modules/checkAVXpe.cpp rename to CMakeModules/checkAVXpe.cpp diff --git a/Modules/checkSSEpe.cpp b/CMakeModules/checkSSEpe.cpp similarity index 100% rename from Modules/checkSSEpe.cpp rename to CMakeModules/checkSSEpe.cpp diff --git a/Modules/compileTestAvx.cpp b/CMakeModules/compileTestAvx.cpp similarity index 100% rename from Modules/compileTestAvx.cpp rename to CMakeModules/compileTestAvx.cpp diff --git a/Modules/compileTestSse.cpp b/CMakeModules/compileTestSse.cpp similarity index 100% rename from Modules/compileTestSse.cpp rename to CMakeModules/compileTestSse.cpp diff --git a/Modules/gcovr b/CMakeModules/gcovr similarity index 100% rename from Modules/gcovr rename to CMakeModules/gcovr diff --git a/Modules/getCpuInfos.cpp b/CMakeModules/getCpuInfos.cpp similarity index 100% rename from Modules/getCpuInfos.cpp rename to CMakeModules/getCpuInfos.cpp diff --git a/README.txt b/README.txt index a364920176e006b6a7c78dac6aef20877fe84c32..38365b3df2c4d8a0694f90d090f80d9b84a551c0 100644 --- a/README.txt +++ b/README.txt @@ -1,8 +1,11 @@ ScalFmm, Inria, Please read the licence. -To compile: --------------------------------------------------- +--------------------------------------------------- + +To compile: +========== # Go to cd scalfmm/Build # Use cmake first by @@ -19,8 +22,11 @@ ccmake .. make # And access executables in scalfmm/Build/Tests/{Release,Debug}/..... -Build the doc --------------------------------------------------- +--------------------------------------------------- + +Build the doc: +============= In scalfmm/Doc you can find several pdf and .tex file about the implementation, kernels and data structure. @@ -32,3 +38,28 @@ make doc browser scalfmm/Build/Doc/html/index.html +--------------------------------------------------- +--------------------------------------------------- + +Getting help and having news from us: +==================================== + +You can subscribe to the scalfmm users mailing list ( scalfmm-public-users@lists.gforge.inria.fr, http://lists.gforge.inria.fr/cgi-bin/mailman/listinfo/scalfmm-public-users ). Very low trafic (~ 2 mails per year) just to know when a news version or an improvement is available. + +Contact the developers at : scalfmm-public-support@lists.gforge.inria.fr + + +--------------------------------------------------- +--------------------------------------------------- + +What inside : +============= +× Src : The Core of Scalfmm is under the Src directory. Users should not need to modify the source. +One can want to implement its own kernel or even its own parallelization whithout modifying the sources. +× Data : example of particles distributions +× Examples : examples of very common usage of Scalfmm +× Doc : should contains the generated Doc +× UTests : contains some unit tests (it can be a good example to understand some features) +× Tests : examples to know how to use scalfmm/put particles in the tree/iterate on the tree... +× Utils : some scripts to work with the data files. + diff --git a/Src/Adaptative/FAbstractAdaptiveKernel.hpp b/Src/Adaptive/FAbstractAdaptiveKernel.hpp similarity index 90% rename from Src/Adaptative/FAbstractAdaptiveKernel.hpp rename to Src/Adaptive/FAbstractAdaptiveKernel.hpp index 25524c9a24de018e379f6ba764df5188d14ae870..57c8679fc26fc9a2df17eb217f102d8fde9b9307 100644 --- a/Src/Adaptative/FAbstractAdaptiveKernel.hpp +++ b/Src/Adaptive/FAbstractAdaptiveKernel.hpp @@ -1,8 +1,8 @@ -#ifndef FABSTRACTADAPTATIVEKERNEL -#define FABSTRACTADAPTATIVEKERNEL +#ifndef FABSTRACTADAPTIVEKERNEL +#define FABSTRACTADAPTIVEKERNEL /** - * This class represent the method that an adaptative kernel must implement. + * This class represent the method that an adaptive kernel must implement. * There are two kinds of operators, the first one represent computation and the others * should return the cretiria to know when the P2M should be performed. */ diff --git a/Src/AdaptiveTree/FAdaptChebSymKernel.hpp b/Src/Adaptive/FAdaptChebSymKernel.hpp similarity index 99% rename from Src/AdaptiveTree/FAdaptChebSymKernel.hpp rename to Src/Adaptive/FAdaptChebSymKernel.hpp index 3e465abc1d56f03f1051edd45191271fc93d9bb5..6cfd793ac75cecb819f31bc2f608644c5fc4385b 100755 --- a/Src/AdaptiveTree/FAdaptChebSymKernel.hpp +++ b/Src/Adaptive/FAdaptChebSymKernel.hpp @@ -19,9 +19,9 @@ #include "Utils/FPoint.hpp" -#include "Adaptative/FAdaptiveCell.hpp" -#include "Adaptative/FAdaptiveKernelWrapper.hpp" -#include "Adaptative/FAbstractAdaptiveKernel.hpp" +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" #include "Kernels/Chebyshev/FChebSymKernel.hpp" class FTreeCoordinate; diff --git a/Src/AdaptiveTree/FAdaptUnifKernel.hpp b/Src/Adaptive/FAdaptUnifKernel.hpp similarity index 99% rename from Src/AdaptiveTree/FAdaptUnifKernel.hpp rename to Src/Adaptive/FAdaptUnifKernel.hpp index 5bd70ae319a4aeadaf25dc6a0a16c50a8ec58e8b..456ed29ef5be299ee38ec7af76707e013029044f 100755 --- a/Src/AdaptiveTree/FAdaptUnifKernel.hpp +++ b/Src/Adaptive/FAdaptUnifKernel.hpp @@ -19,9 +19,9 @@ #include "Utils/FPoint.hpp" -#include "Adaptative/FAdaptiveCell.hpp" -#include "Adaptative/FAdaptiveKernelWrapper.hpp" -#include "Adaptative/FAbstractAdaptiveKernel.hpp" +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" #include "Kernels/Uniform/FUnifKernel.hpp" #include "Kernels/Uniform/FUnifM2LHandler.hpp" diff --git a/Src/Adaptative/FAdaptiveCell.hpp b/Src/Adaptive/FAdaptiveCell.hpp similarity index 100% rename from Src/Adaptative/FAdaptiveCell.hpp rename to Src/Adaptive/FAdaptiveCell.hpp diff --git a/Src/Adaptative/FAdaptiveKernelWrapper.hpp b/Src/Adaptive/FAdaptiveKernelWrapper.hpp similarity index 100% rename from Src/Adaptative/FAdaptiveKernelWrapper.hpp rename to Src/Adaptive/FAdaptiveKernelWrapper.hpp diff --git a/ToDO b/TODO.txt similarity index 100% rename from ToDO rename to TODO.txt diff --git a/Tests/noDist/testAdaptiveChebSymFMM.cpp b/Tests/noDist/testAdaptiveChebSymFMM.cpp index 9e698ebc0845d1e8589dbd8a1cabd1139d9988ff..aef69e823901fd5a2f3f3d13ce8a54e6bcb2aaa6 100644 --- a/Tests/noDist/testAdaptiveChebSymFMM.cpp +++ b/Tests/noDist/testAdaptiveChebSymFMM.cpp @@ -43,9 +43,9 @@ #include "Kernels/Chebyshev/FChebCell.hpp" #include "AdaptiveTree/FAdaptChebSymKernel.hpp" -#include "Adaptative/FAdaptiveCell.hpp" -#include "Adaptative/FAdaptiveKernelWrapper.hpp" -#include "Adaptative/FAbstractAdaptiveKernel.hpp" +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" // #include "Kernels/Interpolation/FInterpMatrixKernel.hpp" #include "Kernels/Chebyshev/FChebCell.hpp" diff --git a/Tests/noDist/testAdaptiveUnifFMM.cpp b/Tests/noDist/testAdaptiveUnifFMM.cpp index c6c66db5bd70823478cd5fc1617e937e9f48b46c..1f0ec20fd3ea742eaa5a42769d0790c97b2d24a8 100644 --- a/Tests/noDist/testAdaptiveUnifFMM.cpp +++ b/Tests/noDist/testAdaptiveUnifFMM.cpp @@ -40,14 +40,14 @@ #include "Components/FSimpleIndexedLeaf.hpp" #include "Kernels/P2P/FP2PParticleContainerIndexed.hpp" -#include "Adaptative/FAdaptiveCell.hpp" -#include "Adaptative/FAdaptiveKernelWrapper.hpp" -#include "Adaptative/FAbstractAdaptiveKernel.hpp" +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" // #include "Kernels/Interpolation/FInterpMatrixKernel.hpp" #include "Kernels/Uniform/FUnifCell.hpp" -#include "AdaptiveTree/FAdaptUnifKernel.hpp" -#include "AdaptiveTree/FAdaptTools.hpp" +#include "Adaptive/FAdaptUnifKernel.hpp" +#include "Adaptive/FAdaptTools.hpp" // // #include "Core/FFmmAlgorithm.hpp" diff --git a/Tests/noDist/testFmmAdaptiveAlgorithm.cpp b/Tests/noDist/testFmmAdaptiveAlgorithm.cpp index bad836406f9a35d7f0a50b3d210d262017996ae3..4c2022634e8aeb3a882e5241a2d4e51383bd8c7b 100644 --- a/Tests/noDist/testFmmAdaptiveAlgorithm.cpp +++ b/Tests/noDist/testFmmAdaptiveAlgorithm.cpp @@ -41,9 +41,9 @@ #include "../../Src/Files/FRandomLoader.hpp" -#include "../../Src/Adaptative/FAdaptiveCell.hpp" -#include "../../Src/Adaptative/FAdaptiveKernelWrapper.hpp" -#include "../../Src/Adaptative/FAbstractAdaptiveKernel.hpp" +#include "../../Src/Adaptive/FAdaptiveCell.hpp" +#include "../../Src/Adaptive/FAdaptiveKernelWrapper.hpp" +#include "../../Src/Adaptive/FAbstractAdaptiveKernel.hpp" template< class CellClass, class ContainerClass> class FAdaptiveTestKernel : public FTestKernels<CellClass, ContainerClass>, public FAbstractAdaptiveKernel<CellClass, ContainerClass> { diff --git a/Tests/noDist/testFmmAdaptiveStats.cpp b/Tests/noDist/testFmmAdaptiveStats.cpp index 68c5521f7972da611ca9d389c4ac6b86a7466608..0c27c98544133b4a43afe938ce515f98458bfbeb 100644 --- a/Tests/noDist/testFmmAdaptiveStats.cpp +++ b/Tests/noDist/testFmmAdaptiveStats.cpp @@ -42,9 +42,9 @@ #include "../../Src/Files/FRandomLoader.hpp" -#include "../../Src/Adaptative/FAdaptiveCell.hpp" -#include "../../Src/Adaptative/FAdaptiveKernelWrapper.hpp" -#include "../../Src/Adaptative/FAbstractAdaptiveKernel.hpp" +#include "../../Src/Adaptive/FAdaptiveCell.hpp" +#include "../../Src/Adaptive/FAdaptiveKernelWrapper.hpp" +#include "../../Src/Adaptive/FAbstractAdaptiveKernel.hpp" template< class CellClass, class ContainerClass> class FAdaptiveStatsKernel : public FAbstractKernels<CellClass, ContainerClass>, public FAbstractAdaptiveKernel<CellClass, ContainerClass> { diff --git a/Src/AdaptiveTree/FAdaptCell.hpp b/ToRemove/AdaptiveTree/FAdaptCell.hpp similarity index 100% rename from Src/AdaptiveTree/FAdaptCell.hpp rename to ToRemove/AdaptiveTree/FAdaptCell.hpp diff --git a/ToRemove/AdaptiveTree/FAdaptChebSymKernel.hpp b/ToRemove/AdaptiveTree/FAdaptChebSymKernel.hpp new file mode 100755 index 0000000000000000000000000000000000000000..6cfd793ac75cecb819f31bc2f608644c5fc4385b --- /dev/null +++ b/ToRemove/AdaptiveTree/FAdaptChebSymKernel.hpp @@ -0,0 +1,469 @@ +#ifndef FADAPTCHEBSYMKERNEL_HPP +#define FADAPTCHEBSYMKERNEL_HPP +// =================================================================================== +// Copyright ScalFmm 2011 INRIA, +// 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". +// =================================================================================== + +#include "Utils/FGlobal.hpp" + +#include "Utils/FPoint.hpp" + +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" +#include "Kernels/Chebyshev/FChebSymKernel.hpp" + +class FTreeCoordinate; + +// ==== CMAKE ===== +// @FUSE_BLAS +// ================ + +// for verbosity only!!! +//#define COUNT_BLOCKED_INTERACTIONS + +// if timings should be logged +//#define LOG_TIMINGS + +/** + * @author O. Coulaud + * @class FAdaptChebSymKernel + * @brief + * Please read the license + * + * This kernels implement the Chebyshev interpolation based FMM operators + * exploiting the symmetries in the far-field. It implements all interfaces + * (P2P, P2M, M2M, M2L, L2L, L2P) which are required by the FFmmAlgorithm and + * FFmmAlgorithmThread. + * + * @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 FAdaptiveChebSymKernel : FChebSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> +, public FAbstractAdaptiveKernel<CellClass, ContainerClass> { + // + typedef FChebSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass; + enum {order = ORDER, + nnodes = TensorTraits<ORDER>::nnodes}; + + const MatrixKernelClass *const MatrixKernel; + +public: + + using KernelBaseClass::P2M; + using KernelBaseClass::M2M; + using KernelBaseClass::M2L; + using KernelBaseClass::finishedLevelM2L; + using KernelBaseClass::L2L; + using KernelBaseClass::L2P; + using KernelBaseClass::P2P; + using KernelBaseClass::P2PRemote; + // /** + // * 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). + // */ + FAdaptiveChebSymKernel(const int inTreeHeight, const FReal inBoxWidth, + const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), MatrixKernel(inMatrixKernel) + {} + // /** Copy constructor */ + FAdaptiveChebSymKernel(const FAdaptiveChebSymKernel& other) + : KernelBaseClass(other), MatrixKernel(other.MatrixKernel) + { } + + // + // /** Destructor */ + ~FAdaptiveChebSymKernel() + { + //this->~KernelBaseClass() ; + } + void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override { + + const FPoint CellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),cellLevel)); + const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + // 1) apply Sy + KernelBaseClass::Interpolator->applyP2M(CellCenter, BoxWidth, + pole->getMultipole(idxRhs), particles); + } + + } + + void M2M(CellClass* const pole, const int poleLevel, const CellClass* const subCell, const int subCellLevel) override { + + const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel)); + const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel))); + + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + + //////////////////////////////////////////////////////////////////////////// + /// p^6 version + // allocate memory + FReal* subChildParentInterpolator = new FReal [nnodes * nnodes]; + + // set child info + FPoint ChildRoots[nnodes], localChildRoots[nnodes]; + FChebTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots); + + // map global position of roots to local position in parent cell + const map_glob_loc map(poleCellCenter, poleCellWidth); + for (unsigned int n=0; n<nnodes; ++n) + map(ChildRoots[n], localChildRoots[n]); + + // assemble child - parent - interpolator + KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // 1) apply Sy (using tensor product M2M with an interpolator computed on the fly) + + // Do NOT reset multipole expansion ! + //FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs)); + + /// p^6 version + FBlas::gemtva(nnodes, nnodes, FReal(1.), + subChildParentInterpolator, + const_cast<FReal*>(subCell->getMultipole(idxRhs)), pole->getMultipole(idxRhs)); + + + }// NVALS + + } + + void P2L(CellClass* const local, const int localLevel, const ContainerClass* const particles) override { + + // Target cell: local + const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel))); + const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + + // interpolation points of target (X) cell + FPoint X[nnodes]; + FChebTensor<order>::setRoots(localCellCenter, localCellWidth, X); + + // read positions + const FReal*const positionsX = particles->getPositions()[0]; + const FReal*const positionsY = particles->getPositions()[1]; + const FReal*const positionsZ = particles->getPositions()[2]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // read physicalValue + const FReal*const physicalValues = particles->getPhysicalValues(); + + // apply P2L + for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){ + + const FPoint y = FPoint(positionsX[idxPart], + positionsY[idxPart], + positionsZ[idxPart]); + + for (unsigned int m=0; m<nnodes; ++m) + local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], y) * physicalValues[idxPart]; + + } + + }// NVALS + + } + + void M2L(CellClass* const local, const int localLevel, const CellClass* const pole, const int poleLevel) override { + + // Source cell: pole + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + + // Target cell: local + const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel))); + const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + + // interpolation points of source (Y) and target (X) cell + FPoint X[nnodes], Y[nnodes]; + FChebTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y); + FChebTensor<order>::setRoots(localCellCenter, localCellWidth, X); + + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // Dense M2L + const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs); + + for (unsigned int m=0; m<nnodes; ++m) + for (unsigned int n=0; n<nnodes; ++n){ + local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n]; + + } + } + } + + void M2P(const CellClass* const pole, const int poleLevel, ContainerClass* const particles) override { + + // Source cell: pole + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + + // interpolation points of source (Y) cell + FPoint Y[nnodes]; + FChebTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y); + + // read positions + const FReal*const positionsX = particles->getPositions()[0]; + const FReal*const positionsY = particles->getPositions()[1]; + const FReal*const positionsZ = particles->getPositions()[2]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs); + + // apply M2P + for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){ + + const FPoint x = FPoint(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]); + + FReal targetValue=0.; + { + for (unsigned int n=0; n<nnodes; ++n) + targetValue += + MatrixKernel->evaluate(x, Y[n]) * MultipoleExpansion[/*idxLhs*nnodes+*/n]; + } + + // get potential + FReal*const potentials = particles->getPotentials(/*idxPot*/); + // add contribution to potential + potentials[idxPart] += (targetValue); + + }// Particles + + }// NVALS + + } + + void L2L(const CellClass* const local, const int localLevel, CellClass* const subCell, const int subCellLevel) override { + + const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel)); + const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel))); + + const FPoint localCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + const FReal localWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,localLevel))); + + //////////////////////////////////////////////////////////////////////////// + /// p^6 version + // allocate memory + FReal* subChildParentInterpolator = new FReal [nnodes * nnodes]; + + // set child info + FPoint ChildRoots[nnodes], localChildRoots[nnodes]; + FChebTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots); + + // map global position of roots to local position in parent cell + const map_glob_loc map(localCenter, localWidth); + for (unsigned int n=0; n<nnodes; ++n) + map(ChildRoots[n], localChildRoots[n]); + + // assemble child - parent - interpolator + KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // 2) apply Sx + + /// p^6 version + FBlas::gemva(nnodes, nnodes, FReal(1.), + subChildParentInterpolator, + const_cast<FReal*>(local->getLocal(idxRhs)), subCell->getLocal(idxRhs)); + + }// NVALS + + } + + void L2P(const CellClass* const local, const int cellLevel, ContainerClass* const particles) override { + + const FPoint CellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),cellLevel)); + const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // 2.a) apply Sx + KernelBaseClass::Interpolator->applyL2P(CellCenter, BoxWidth, + local->getLocal(idxRhs), particles); + + // 2.b) apply Px (grad Sx) + KernelBaseClass::Interpolator->applyL2PGradient(CellCenter, BoxWidth, + local->getLocal(idxRhs), particles); + + } + + } + + void P2P(ContainerClass* target, const ContainerClass* sources) override { +// long long int*const particlesAttributes = target->getDataDown(); +// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){ +// particlesAttributes[idxPart] += sources->getNbParticles(); +// } + } + + bool preferP2M(const ContainerClass* const particles) override { + return particles->getNbParticles() < 10; + } + bool preferP2M(const int /*atLevel*/, const ContainerClass*const particles[], const int nbContainers) override { + int counterParticles = 0; + for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){ + counterParticles += particles[idxContainer]->getNbParticles(); + } + return counterParticles < 10; + } +}; + +// +//template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1> +//class FAdaptChebSymKernel +// : public FChebSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> +//{ +// typedef FChebSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass; +// +//#ifdef LOG_TIMINGS +// FTic time; +// FReal t_m2l_1, t_m2l_2, t_m2l_3; +//#endif +// +//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). +// */ +// FAdaptChebSymKernel(const int inTreeHeight, +// const FReal inBoxWidth, +// const FPoint& inBoxCenter) +//: KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter) +//{ +// +//#ifdef LOG_TIMINGS +// t_m2l_1 = FReal(0.); +// t_m2l_2 = FReal(0.); +// t_m2l_3 = FReal(0.); +//#endif +//} +// +// +// /** Copy constructor */ +// FAdaptChebSymKernel(const FAdaptChebSymKernel& other) +// : KernelBaseClass(other) +// { } +// +// +// +// /** Destructor */ +// ~FAdaptChebSymKernel() +// { +// this->~KernelBaseClass() ; +//#ifdef LOG_TIMINGS +// std::cout << "- Permutation took " << t_m2l_1 << "s" +// << "\n- GEMMT and GEMM took " << t_m2l_2 << "s" +// << "\n- Unpermutation took " << t_m2l_3 << "s" +// << std::endl; +//#endif +// } +// +// +// void P2MAdapt(CellClass* const ParentCell, const int &level) +// { +// const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(ParentCell->getCoordinate())); +// const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-level); +// // +// for(int i = 0 ; i <ParentCell->getLeavesSize(); ++i ){ +// // +// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth, +// ParentCell->getMultipole(idxRhs), ParentCell->getLeaf(i)->getSrc()); +// } +// } +// } +// void M2MAdapt(CellClass* const FRestrict ParentCell, const int &TreeLevel, const int &numberOfM2M, +// const int * FRestrict ChildLevel , const CellClass*const FRestrict *const FRestrict ChildCells) +// { +// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// // // apply Sy +// for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ +// if (ChildCells[ChildIndex]){ +// // KernelBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs)); +// } +// } +// } +// } +// +// +// +// void M2L(CellClass* const FRestrict TargetCell, +// const CellClass* SourceCells[343], +// const int /*NumSourceCells*/, +// const int TreeLevel) +// { +// +// } +// +// +// void L2L(const CellClass* const FRestrict ParentCell, +// CellClass* FRestrict *const FRestrict ChildCells, +// const int /*TreeLevel*/) +// { +// // for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// // // apply Sx +// // for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ +// // if (ChildCells[ChildIndex]){ +// // AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs)); +// // } +// // } +// // } +// } +// +// void L2P(const CellClass* const LeafCell, +// ContainerClass* const TargetParticles) +// { +// KernelBaseClass::L2P(LeafCell,TargetParticles) ; +// } +// +// // 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 //FADAPTCHEBSYMKERNELS_HPP + +// [--END--] diff --git a/Src/AdaptiveTree/FAdaptSeqAlgorithm.hpp b/ToRemove/AdaptiveTree/FAdaptSeqAlgorithm.hpp similarity index 100% rename from Src/AdaptiveTree/FAdaptSeqAlgorithm.hpp rename to ToRemove/AdaptiveTree/FAdaptSeqAlgorithm.hpp diff --git a/Src/AdaptiveTree/FAdaptTools.hpp b/ToRemove/AdaptiveTree/FAdaptTools.hpp similarity index 100% rename from Src/AdaptiveTree/FAdaptTools.hpp rename to ToRemove/AdaptiveTree/FAdaptTools.hpp diff --git a/ToRemove/AdaptiveTree/FAdaptUnifKernel.hpp b/ToRemove/AdaptiveTree/FAdaptUnifKernel.hpp new file mode 100755 index 0000000000000000000000000000000000000000..456ed29ef5be299ee38ec7af76707e013029044f --- /dev/null +++ b/ToRemove/AdaptiveTree/FAdaptUnifKernel.hpp @@ -0,0 +1,593 @@ +#ifndef FADAPTUNIFKERNEL_HPP +#define FADAPTUNIFKERNEL_HPP +// =================================================================================== +// Copyright ScalFmm 2011 INRIA, +// 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". +// =================================================================================== + +#include "Utils/FGlobal.hpp" + +#include "Utils/FPoint.hpp" + +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" +#include "Kernels/Uniform/FUnifKernel.hpp" +#include "Kernels/Uniform/FUnifM2LHandler.hpp" + +class FTreeCoordinate; + +// ==== CMAKE ===== +// @FUSE_BLAS +// ================ + +// for verbosity only!!! +//#define COUNT_BLOCKED_INTERACTIONS + +// if timings should be logged +//#define LOG_TIMINGS + +/** + * @author O. Coulaud + * @class FAdaptUnifKernel + * @brief + * Please read the license + * + * This kernels implement the Lagrange interpolation based FMM operators. + * It implements all interfaces (P2P, P2M, M2M, M2L, L2L, L2P) which are + * required by the FFmmAlgorithm and FFmmAlgorithmThread. + * + * @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 FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> +, public FAbstractAdaptiveKernel<CellClass, ContainerClass> { + // + typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass; + + enum {order = ORDER, + nnodes = TensorTraits<ORDER>::nnodes}; + + /// Needed for M2L operator +// // If we choose to pre-assemble adaptive M2L operators +// // then we need to provide an adaptive M2L handler +// // and the transfer in Fourier space is straightforward. +// typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass; +// const M2LHandlerClass M2LHandler; + + const MatrixKernelClass *const MatrixKernel; + +public: + + using KernelBaseClass::P2M; + using KernelBaseClass::M2M; + using KernelBaseClass::M2L; + using KernelBaseClass::finishedLevelM2L; + using KernelBaseClass::L2L; + using KernelBaseClass::L2P; + using KernelBaseClass::P2P; + using KernelBaseClass::P2PRemote; + // /** + // * 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). + // */ + FAdaptiveUnifKernel(const int inTreeHeight, const FReal inBoxWidth, + const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel)/*, M2LHandler(inMatrixKernel, inTreeHeight, inBoxWidth)*/, MatrixKernel(inMatrixKernel) + {} + // /** Copy constructor */ + FAdaptiveUnifKernel(const FAdaptiveUnifKernel& other) + : KernelBaseClass(other)/*, M2LHandler(other.M2LHandler)*/, MatrixKernel(other.MatrixKernel) + { } + + // + // /** Destructor */ + ~FAdaptiveUnifKernel() + { + //this->~KernelBaseClass() ; + } + void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override { + + const FPoint CellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),cellLevel)); + const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + // 1) apply Sy + KernelBaseClass::Interpolator->applyP2M(CellCenter, BoxWidth, + pole->getMultipole(idxRhs), particles); +// // 2) apply Discrete Fourier Transform +// M2LHandler.applyZeroPaddingAndDFT(pole->getMultipole(idxRhs), +// pole->getTransformedMultipole(idxRhs)); + + } + + } + + void M2M(CellClass* const pole, const int poleLevel, const CellClass* const subCell, const int subCellLevel) override { + + const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel)); + const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel))); + + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + +// //////////////////////////////////////////////////////////////////////////// +// /// p^6 version +// // allocate memory +// FReal* subChildParentInterpolator = new FReal [nnodes * nnodes]; +// +// // set child info +// FPoint ChildRoots[nnodes], localChildRoots[nnodes]; +// FUnifTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots); +// +// // map global position of roots to local position in parent cell +// const map_glob_loc map(poleCellCenter, poleCellWidth); +// for (unsigned int n=0; n<nnodes; ++n) +// map(ChildRoots[n], localChildRoots[n]); +// +// // assemble child - parent - interpolator +// KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator); + + + //////////////////////////////////////////////////////////////////////////// + /// p^4 version + + // Set sub-child coords + FReal globalChildCoords[3][ORDER]; + FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords); + + // Map global position of sub-child nodes to [-1,1] + FReal localChildCoords[3][ORDER]; + const map_glob_loc map(poleCellCenter, poleCellWidth); + FPoint localChildPoints; + for (unsigned int n=0; n<ORDER; ++n) { + map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints); + localChildCoords[0][n] = localChildPoints.getX(); + localChildCoords[1][n] = localChildPoints.getY(); + localChildCoords[2][n] = localChildPoints.getZ(); + } + + // assemble interpolator + FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER]; + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator); + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER); + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER); + + // get permutation operators + unsigned int perm[3][nnodes]; + for (unsigned int i=0;i<3; ++i) + for (unsigned int n=0; n<nnodes; ++n) + perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // 1) apply Sy (using tensor product M2M with an interpolator computed on the fly) + + // Do NOT reset multipole expansion ! + //FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs)); + +// /// p^6 version +// FBlas::gemtva(nnodes, nnodes, FReal(1.), +// subChildParentInterpolator, +// const_cast<FReal*>(subCell->getMultipole(idxRhs)), pole->getMultipole(idxRhs)); + + /// p^4 version + FReal Exp[nnodes], PermExp[nnodes]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator, ORDER, + const_cast<FReal*>(subCell->getMultipole(idxRhs)), ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator + 2 * ORDER*ORDER, ORDER, + Exp, ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator + 1 * ORDER*ORDER, ORDER, + Exp, ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) pole->getMultipole(idxRhs)[perm[2][n]] += PermExp[n]; + + + +// // 2) Apply Discete Fourier Transform +// M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs), +// ParentCell->getTransformedMultipole(idxRhs)); + } + } + + void P2L(CellClass* const local, const int localLevel, const ContainerClass* const particles) override { + + // Target cell: local + const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel))); + const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + + // interpolation points of target (X) cell + FPoint X[nnodes]; + FUnifTensor<order>::setRoots(localCellCenter, localCellWidth, X); + + // read positions + const FReal*const positionsX = particles->getPositions()[0]; + const FReal*const positionsY = particles->getPositions()[1]; + const FReal*const positionsZ = particles->getPositions()[2]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // read physicalValue + const FReal*const physicalValues = particles->getPhysicalValues(); + + // apply P2L + for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){ + + const FPoint y = FPoint(positionsX[idxPart], + positionsY[idxPart], + positionsZ[idxPart]); + + for (unsigned int m=0; m<nnodes; ++m) + local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], y) * physicalValues[idxPart]; + + } + + }// NVALS + + } + + void M2L(CellClass* const local, const int localLevel, const CellClass* const pole, const int poleLevel) override { + + // Source cell: pole + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + + // Target cell: local + const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel))); + const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + + // interpolation points of source (Y) and target (X) cell + FPoint X[nnodes], Y[nnodes]; + FUnifTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y); + FUnifTensor<order>::setRoots(localCellCenter, localCellWidth, X); + + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + // Dense M2L + const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs); + + for (unsigned int m=0; m<nnodes; ++m) + for (unsigned int n=0; n<nnodes; ++n){ + local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n]; + + } + } + } + + void M2P(const CellClass* const pole, const int poleLevel, ContainerClass* const particles) override { + + // Source cell: pole + const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel))); + const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel)); + + // interpolation points of source (Y) cell + FPoint Y[nnodes]; + FUnifTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y); + + // read positions + const FReal*const positionsX = particles->getPositions()[0]; + const FReal*const positionsY = particles->getPositions()[1]; + const FReal*const positionsZ = particles->getPositions()[2]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + + const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs); + + // apply M2P + for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){ + + const FPoint x = FPoint(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]); + + FReal targetValue=0.; + { + for (unsigned int n=0; n<nnodes; ++n) + targetValue += + MatrixKernel->evaluate(x, Y[n]) * MultipoleExpansion[/*idxLhs*nnodes+*/n]; + } + + // get potential + FReal*const potentials = particles->getPotentials(/*idxPot*/); + // add contribution to potential + potentials[idxPart] += (targetValue); + + }// Particles + + }// NVALS + + } + + void L2L(const CellClass* const local, const int localLevel, CellClass* const subCell, const int subCellLevel) override { + + const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel)); + const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel))); + + const FPoint localCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel)); + const FReal localWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,localLevel))); + +// //////////////////////////////////////////////////////////////////////////// +// /// p^6 version +// // allocate memory +// FReal* subChildParentInterpolator = new FReal [nnodes * nnodes]; +// +// // set child info +// FPoint ChildRoots[nnodes], localChildRoots[nnodes]; +// FUnifTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots); +// +// // map global position of roots to local position in parent cell +// const map_glob_loc map(localCenter, localWidth); +// for (unsigned int n=0; n<nnodes; ++n) +// map(ChildRoots[n], localChildRoots[n]); +// +// // assemble child - parent - interpolator +// KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator); + + //////////////////////////////////////////////////////////////////////////// + /// p^4 version + // Set sub-child coords + FReal globalChildCoords[3][ORDER]; + FUnifTensor<order>::setPolynomialsRoots(subCellCenter, subCellWidth, globalChildCoords); + + // Map global position of sub-child nodes to [-1,1] + FReal localChildCoords[3][ORDER]; + const map_glob_loc map(localCenter, localWidth); + FPoint localChildPoints; + for (unsigned int n=0; n<ORDER; ++n) { + map(FPoint(globalChildCoords[0][n],globalChildCoords[1][n],globalChildCoords[2][n]), localChildPoints); + localChildCoords[0][n] = localChildPoints.getX(); + localChildCoords[1][n] = localChildPoints.getY(); + localChildCoords[2][n] = localChildPoints.getZ(); + } + + // assemble interpolator + FReal* subChildParentInterpolator = new FReal [3 * ORDER*ORDER]; + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[0], subChildParentInterpolator); + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[1], subChildParentInterpolator + 1 * ORDER*ORDER); + KernelBaseClass::Interpolator->assembleInterpolator(ORDER, localChildCoords[2], subChildParentInterpolator + 2 * ORDER*ORDER); + + // get permutation operators + unsigned int perm[3][nnodes]; + for (unsigned int i=0;i<3; ++i) + for (unsigned int n=0; n<nnodes; ++n) + perm[i][n] = KernelBaseClass::Interpolator->getPermutationsM2ML2L(i)[n]; + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// // 1) Apply Inverse Discete Fourier Transform +// M2LHandler.unapplyZeroPaddingAndDFT(local->getTransformedLocal(idxRhs), +// const_cast<CellClass*>(local)->getLocal(idxRhs)); + + // 2) apply Sx + +// /// p^6 version +// FBlas::gemva(nnodes, nnodes, FReal(1.), +// subChildParentInterpolator, +// const_cast<FReal*>(local->getLocal(idxRhs)), subCell->getLocal(idxRhs)); + + /// p^4 version + FReal Exp[nnodes], PermExp[nnodes]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator, ORDER, + const_cast<FReal*>(local->getLocal(idxRhs)), ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator + 2 * ORDER*ORDER, ORDER, + Exp, ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]]; + // ORDER*ORDER*ORDER * (2*ORDER-1) + FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), + subChildParentInterpolator + 1 * ORDER*ORDER, ORDER, + Exp, ORDER, PermExp, ORDER); + + for (unsigned int n=0; n<nnodes; ++n) subCell->getLocal(idxRhs)[perm[2][n]] += PermExp[n]; + // total flops count: 3 * ORDER*ORDER*ORDER * (2*ORDER-1) + + } + + } + + void L2P(const CellClass* const local, const int cellLevel, ContainerClass* const particles) override { + + const FPoint CellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),cellLevel)); + const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel); + + for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ + +// // 1) Apply Inverse Discete Fourier Transform +// M2LHandler.unapplyZeroPaddingAndDFT(local->getTransformedLocal(idxRhs), +// const_cast<CellClass*>(local)->getLocal(idxRhs)); + + // 2.a) apply Sx + KernelBaseClass::Interpolator->applyL2P(CellCenter, BoxWidth, + local->getLocal(idxRhs), particles); + + // 2.b) apply Px (grad Sx) + KernelBaseClass::Interpolator->applyL2PGradient(CellCenter, BoxWidth, + local->getLocal(idxRhs), particles); + + } + } + + void P2P(ContainerClass* target, const ContainerClass* sources) override { +// long long int*const particlesAttributes = target->getDataDown(); +// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){ +// particlesAttributes[idxPart] += sources->getNbParticles(); +// } + } + + bool preferP2M(const ContainerClass* const particles) override { + return particles->getNbParticles() < 10; + } + bool preferP2M(const int /*atLevel*/, const ContainerClass*const particles[], const int nbContainers) override { + int counterParticles = 0; + for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){ + counterParticles += particles[idxContainer]->getNbParticles(); + } + return counterParticles < 10; + } +}; + +// +//template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1> +//class FAdaptUnifKernel +// : public FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> +//{ +// typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass; +// +//#ifdef LOG_TIMINGS +// FTic time; +// FReal t_m2l_1, t_m2l_2, t_m2l_3; +//#endif +// +//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). +// */ +// FAdaptUnifKernel(const int inTreeHeight, +// const FReal inBoxWidth, +// const FPoint& inBoxCenter) +//: KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter) +//{ +// +//#ifdef LOG_TIMINGS +// t_m2l_1 = FReal(0.); +// t_m2l_2 = FReal(0.); +// t_m2l_3 = FReal(0.); +//#endif +//} +// +// +// /** Copy constructor */ +// FAdaptUnifKernel(const FAdaptUnifKernel& other) +// : KernelBaseClass(other) +// { } +// +// +// +// /** Destructor */ +// ~FAdaptUnifKernel() +// { +// this->~KernelBaseClass() ; +//#ifdef LOG_TIMINGS +// std::cout << "- Permutation took " << t_m2l_1 << "s" +// << "\n- GEMMT and GEMM took " << t_m2l_2 << "s" +// << "\n- Unpermutation took " << t_m2l_3 << "s" +// << std::endl; +//#endif +// } +// +// +// void P2MAdapt(CellClass* const ParentCell, const int &level) +// { +// const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(ParentCell->getCoordinate())); +// const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-level); +// // +// for(int i = 0 ; i <ParentCell->getLeavesSize(); ++i ){ +// // +// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth, +// ParentCell->getMultipole(idxRhs), ParentCell->getLeaf(i)->getSrc()); +// } +// } +// } +// void M2MAdapt(CellClass* const FRestrict ParentCell, const int &TreeLevel, const int &numberOfM2M, +// const int * FRestrict ChildLevel , const CellClass*const FRestrict *const FRestrict ChildCells) +// { +// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// // // apply Sy +// for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ +// if (ChildCells[ChildIndex]){ +// // KernelBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs)); +// } +// } +// } +// } +// +// +// +// void M2L(CellClass* const FRestrict TargetCell, +// const CellClass* SourceCells[343], +// const int /*NumSourceCells*/, +// const int TreeLevel) +// { +// +// } +// +// +// void L2L(const CellClass* const FRestrict ParentCell, +// CellClass* FRestrict *const FRestrict ChildCells, +// const int /*TreeLevel*/) +// { +// // for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ +// // // apply Sx +// // for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ +// // if (ChildCells[ChildIndex]){ +// // AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs)); +// // } +// // } +// // } +// } +// +// void L2P(const CellClass* const LeafCell, +// ContainerClass* const TargetParticles) +// { +// KernelBaseClass::L2P(LeafCell,TargetParticles) ; +// } +// +// // 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 //FADAPTUNIFKERNELS_HPP + +// [--END--] diff --git a/Tests/Olivier/ChebyshevDirectPeriodic.cpp b/ToRemove/Olivier/ChebyshevDirectPeriodic.cpp similarity index 100% rename from Tests/Olivier/ChebyshevDirectPeriodic.cpp rename to ToRemove/Olivier/ChebyshevDirectPeriodic.cpp diff --git a/Tests/Olivier/DirectEwaldcmp.cpp b/ToRemove/Olivier/DirectEwaldcmp.cpp similarity index 100% rename from Tests/Olivier/DirectEwaldcmp.cpp rename to ToRemove/Olivier/DirectEwaldcmp.cpp diff --git a/Tests/Olivier/LagrangeInterpolation.cpp b/ToRemove/Olivier/LagrangeInterpolation.cpp similarity index 100% rename from Tests/Olivier/LagrangeInterpolation.cpp rename to ToRemove/Olivier/LagrangeInterpolation.cpp diff --git a/Tests/Olivier/adaptiveOctree.cpp b/ToRemove/Olivier/adaptiveOctree.cpp similarity index 100% rename from Tests/Olivier/adaptiveOctree.cpp rename to ToRemove/Olivier/adaptiveOctree.cpp diff --git a/Tests/Olivier/checkCPP.cpp b/ToRemove/Olivier/checkCPP.cpp similarity index 100% rename from Tests/Olivier/checkCPP.cpp rename to ToRemove/Olivier/checkCPP.cpp diff --git a/Tests/Olivier/checkLoader.cpp b/ToRemove/Olivier/checkLoader.cpp similarity index 100% rename from Tests/Olivier/checkLoader.cpp rename to ToRemove/Olivier/checkLoader.cpp diff --git a/Tests/Olivier/statAdapt.cpp b/ToRemove/Olivier/statAdapt.cpp similarity index 100% rename from Tests/Olivier/statAdapt.cpp rename to ToRemove/Olivier/statAdapt.cpp diff --git a/Tests/Olivier/testFMMEwald.cpp b/ToRemove/Olivier/testFMMEwald.cpp similarity index 100% rename from Tests/Olivier/testFMMEwald.cpp rename to ToRemove/Olivier/testFMMEwald.cpp diff --git a/Tests/Olivier/testSphericalKernel.cpp b/ToRemove/Olivier/testSphericalKernel.cpp similarity index 100% rename from Tests/Olivier/testSphericalKernel.cpp rename to ToRemove/Olivier/testSphericalKernel.cpp diff --git a/Tests/TestPara/FChebOrder.cpp b/ToRemove/TestPara/FChebOrder.cpp similarity index 100% rename from Tests/TestPara/FChebOrder.cpp rename to ToRemove/TestPara/FChebOrder.cpp diff --git a/email-sorted b/ToRemove/email-sorted similarity index 100% rename from email-sorted rename to ToRemove/email-sorted diff --git a/utestChebyshevDirectPeriodic.cpp b/ToRemove/utestChebyshevDirectPeriodic.cpp similarity index 100% rename from utestChebyshevDirectPeriodic.cpp rename to ToRemove/utestChebyshevDirectPeriodic.cpp