FUnifDenseKernel.hpp 8.56 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
16 17
// Keep in private GIT

18 19 20 21
#ifndef FUNIFDENSEKERNEL_HPP
#define FUNIFDENSEKERNEL_HPP

#include "../../Utils/FGlobal.hpp"
22

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
#include "../../Utils/FSmartPointer.hpp"

#include "./FAbstractUnifKernel.hpp"
#include "./FUnifM2LHandler.hpp"

class FTreeCoordinate;

/**
 * @author Pierre Blanchard (pierre.blanchard@inria.fr)
 * @class FUnifDenseKernel
 * @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. 
 * This is the dense version of the kernel. The transfer are done in real space
 * and not in Fourier space. 
 *
 * @tparam CellClass Type of cell
 * @tparam ContainerClass Type of container to store particles
 * @tparam MatrixKernelClass Type of matrix kernel function
 * @tparam ORDER Lagrange interpolation order
 */
47
template < class FReal, class CellClass,	class ContainerClass,	class MatrixKernelClass, int ORDER, int NVALS = 1>
48
class FUnifDenseKernel
49
  : public FAbstractUnifKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
50
{
51
    // private types
52
    typedef FUnifM2LHandler<FReal, ORDER,MatrixKernelClass::Type> M2LHandlerClass;
53

54
    // using from
55
    typedef FAbstractUnifKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
56
    AbstractBaseClass;
57

58 59
    /// Needed for P2P and M2L operators
    const MatrixKernelClass *const MatrixKernel;
60

61 62
    /// Needed for M2L operator
    const M2LHandlerClass M2LHandler;
63 64 65


public:
66 67 68 69 70 71 72
    /**
    * 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).
    */
    FUnifDenseKernel(const int inTreeHeight,
                     const FReal inBoxWidth,
73
                     const FPoint<FReal>& inBoxCenter,
74
                     const MatrixKernelClass *const inMatrixKernel)
75
    : FAbstractUnifKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
76 77 78 79
      MatrixKernel(inMatrixKernel),
      M2LHandler(MatrixKernel,
                 inTreeHeight,
                 inBoxWidth) 
80
    { }
81 82


83 84 85
    void P2M(CellClass* const LeafCell,
             const ContainerClass* const SourceParticles)
    {
86
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
87 88 89 90 91 92
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 1) apply Sy
            AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                      LeafCell->getMultipole(idxRhs), SourceParticles);
        }
    }
93

94 95 96 97 98 99 100 101 102 103 104 105 106

    void M2M(CellClass* const FRestrict ParentCell,
             const CellClass*const FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
                if (ChildCells[ChildIndex]){
                    AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
                                                              ParentCell->getMultipole(idxRhs));
                }
            }
        }
107 108
    }

109 110
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int /*neighborPositions*/[], const int inSize, const int TreeLevel)  override {
111
        const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
112

113
        // interpolation points of source (Y) and target (X) cell
114
        FPoint<FReal> X[AbstractBaseClass::nnodes], Y[AbstractBaseClass::nnodes];
115
        FUnifTensor<FReal,ORDER>::setRoots(AbstractBaseClass::getCellCenter(TargetCell->getCoordinate(),TreeLevel), CellWidth, X);
116

117
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
118

119 120
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                FUnifTensor<FReal,ORDER>::setRoots(AbstractBaseClass::getCellCenter(SourceCells[idxExistingNeigh]->getCoordinate(),TreeLevel), CellWidth, Y);
121

122 123 124 125
                for (unsigned int m=0; m<AbstractBaseClass::nnodes; ++m)
                    for (unsigned int n=0; n<AbstractBaseClass::nnodes; ++n){
                        TargetCell->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * SourceCells[idxExistingNeigh]->getMultipole(idxRhs)[n];
                    }
126

127
            }
128 129
        }
    }
130 131 132 133 134 135 136 137 138 139 140 141 142


    void L2L(const CellClass* const FRestrict ParentCell,
             CellClass* FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 2) apply Sx
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
                if (ChildCells[ChildIndex]){
                    AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
                }
            }
143 144 145
        }
    }

146 147 148
    void L2P(const CellClass* const LeafCell,
             ContainerClass* const TargetParticles)
    {
149
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
150

151
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
152

153 154 155
            // 2.a) apply Sx
            AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                      LeafCell->getLocal(idxRhs), TargetParticles);
156

157 158 159 160 161 162 163
            // 2.b) apply Px (grad Sx)
            AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                              LeafCell->getLocal(idxRhs), TargetParticles);

        }
    }

164
    void P2P(const FTreeCoordinate& inPosition,
165
             ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict inSources,
166 167
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
168 169 170 171 172 173 174 175 176
        if(inTargets == inSources){
            P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PInner(inTargets,MatrixKernel);
        }
        else{
            const ContainerClass* const srcPtr[1] = {inSources};
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,srcPtr,1,MatrixKernel);
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
        }
177 178 179 180 181 182
    }

    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
             ContainerClass* const FRestrict inTargets,
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
183 184 185 186 187 188
        int nbNeighborsToCompute = 0;
        while(nbNeighborsToCompute < inSize
              && neighborPositions[nbNeighborsToCompute] < 14){
            nbNeighborsToCompute += 1;
        }
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2P(inTargets,inNeighbors,nbNeighborsToCompute,MatrixKernel);
189
    }
190 191


192 193
    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
194
                   const ContainerClass* const inNeighbors[], const int /*neighborPositions*/[],
195 196
                   const int inSize) override {
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
197 198 199 200 201 202 203 204
    }

};


#endif //FUNIFDENSEKERNEL_HPP

// [--END--]