FUnifKernel.hpp 10.6 KB
Newer Older
1
// ===================================================================================
COULAUD Olivier's avatar
COULAUD Olivier committed
2
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
3 4 5 6 7 8 9 10 11 12 13 14 15
// 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 18
// Keep in private GIT
// @SCALFMM_PRIVATE

19 20 21
#ifndef FUNIFKERNEL_HPP
#define FUNIFKERNEL_HPP

COULAUD Olivier's avatar
COULAUD Olivier committed
22
#include "Utils/FGlobal.hpp"
23

COULAUD Olivier's avatar
COULAUD Olivier committed
24
#include "Utils/FSmartPointer.hpp"
25

COULAUD Olivier's avatar
COULAUD Olivier committed
26 27
#include "FAbstractUnifKernel.hpp"
#include "FUnifM2LHandler.hpp"
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

class FTreeCoordinate;

/**
 * @author Pierre Blanchard (pierre.blanchard@inria.fr)
 * @class FUnifKernel
 * @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 Lagrange interpolation order
 */
46
template < class FReal, class CellClass, class ContainerClass,   class MatrixKernelClass, int ORDER, int NVALS = 1>
47
class FUnifKernel
48
  : public FAbstractUnifKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
49
{
50
    // private types
51
    typedef FUnifM2LHandler<FReal, ORDER,MatrixKernelClass::Type> M2LHandlerClass;
52

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

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

60 61
    /// Needed for M2L operator
    const M2LHandlerClass M2LHandler;
62

63 64
    /// Leaf level separation criterion
    const int LeafLevelSeparationCriterion;
65 66

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


87 88 89
    void P2M(CellClass* const LeafCell,
             const ContainerClass* const SourceParticles)
    {
90
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
91 92 93
        // 1) apply Sy
        AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                  LeafCell->getMultipole(0), SourceParticles);
94

95
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
96

97 98 99
            // 2) apply Discrete Fourier Transform
            M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxRhs), 
                                              LeafCell->getTransformedMultipole(idxRhs));
100

101
        }
102
    }
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120


    void M2M(CellClass* const FRestrict ParentCell,
             const CellClass*const FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 1) apply Sy
            //FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
                if (ChildCells[ChildIndex]){
                    AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
                                                              ParentCell->getMultipole(idxRhs));
                }
            }
            // 2) Apply Discete Fourier Transform
            M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs), 
                                              ParentCell->getTransformedMultipole(idxRhs));
121 122
        }
    }
123 124


125 126
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int neighborPositions[], const int inSize, const int TreeLevel)  override {
127 128 129 130
        const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
        const FReal scale(MatrixKernel->getScaleFactor(CellWidth));

        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
131
            FComplex<FReal> *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);
132

133 134 135 136 137
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idxNeigh = neighborPositions[idxExistingNeigh];
                M2LHandler.applyFC(idxNeigh, TreeLevel, scale,
                                   SourceCells[idxExistingNeigh]->getTransformedMultipole(idxRhs),
                                   TransformedLocalExpansion);
138
            }
139 140
        }
    }
141 142 143 144 145 146 147


    void L2L(const CellClass* const FRestrict ParentCell,
             CellClass* FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
148

149
            // 1) Apply Inverse Discete Fourier Transform
150
            FReal localExp[AbstractBaseClass::nnodes];
151
            M2LHandler.unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxRhs),
152 153 154
                                                localExp);
            FBlas::add(AbstractBaseClass::nnodes,const_cast<FReal*>(ParentCell->getLocal(idxRhs)),localExp);

155 156 157
            // 2) apply Sx
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
                if (ChildCells[ChildIndex]){
158
                    AbstractBaseClass::Interpolator->applyL2L(ChildIndex, localExp, ChildCells[ChildIndex]->getLocal(idxRhs));
159 160
                }
            }
161 162 163
        }
    }

164 165 166
    void L2P(const CellClass* const LeafCell,
             ContainerClass* const TargetParticles)
    {
167
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
168

169 170
        FReal localExp[NVALS*AbstractBaseClass::nnodes];

171
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
172

173 174
            // 1)  Apply Inverse Discete Fourier Transform
            M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs), 
175 176
                                                localExp + idxRhs*AbstractBaseClass::nnodes);
            FBlas::add(AbstractBaseClass::nnodes,const_cast<FReal*>(LeafCell->getLocal(idxRhs)),localExp + idxRhs*AbstractBaseClass::nnodes);
177

178
        }
179

180 181
        // 2.a) apply Sx
        AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
182
                                                  localExp, TargetParticles);
183

184 185
        // 2.b) apply Px (grad Sx)
        AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
186
                                                          localExp, TargetParticles);
187 188


189
    }
190

191
    void P2P(const FTreeCoordinate& inPosition,
192
             ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict inSources,
193 194
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
195
        // Standard FMM separation criterion, i.e. max 27 neighbor clusters per leaf
196 197 198 199 200 201 202 203 204 205
        if(LeafLevelSeparationCriterion==1) {
            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);
            }
206
        }
207
        // Nearfield interactions are only computed within the target leaf
208 209 210
        else if(LeafLevelSeparationCriterion==0){
            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
        }
211
        // If criterion equals -1 then no P2P need to be performed.
212
    }
213

214 215 216 217 218 219 220 221 222 223 224
    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
             ContainerClass* const FRestrict inTargets,
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
        int nbNeighborsToCompute = 0;
        while(nbNeighborsToCompute < inSize
              && neighborPositions[nbNeighborsToCompute] < 14){
            nbNeighborsToCompute += 1;
        }
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2P(inTargets,inNeighbors,nbNeighborsToCompute,MatrixKernel);
    }
225

226 227
    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
228
                   const ContainerClass* const inNeighbors[], const int /*neighborPositions*/[],
229
                   const int inSize) override {
230 231
        // Standard FMM separation criterion, i.e. max 27 neighbor clusters per leaf
        if(LeafLevelSeparationCriterion==1) 
232
            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
233 234 235 236
        // Nearfield interactions are only computed within the target leaf
        if(LeafLevelSeparationCriterion==0) 
            DirectInteractionComputer<FReal,MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,0,MatrixKernel);
        // If criterion equals -1 then no P2P need to be performed.
237
    }
238 239 240 241 242 243 244

};


#endif //FUNIFKERNEL_HPP

// [--END--]