FChebKernel.hpp 11.5 KB
Newer Older
BRAMAS Berenger's avatar
BRAMAS Berenger committed
1
// ===================================================================================
2
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
3 4 5 6 7 8 9 10 11 12 13 14
// 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".
BRAMAS Berenger's avatar
BRAMAS Berenger committed
15
// ===================================================================================
16 17 18
#ifndef FCHEBKERNEL_HPP
#define FCHEBKERNEL_HPP

19
#include "../../Utils/FGlobal.hpp"
20

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

#include "./FAbstractChebKernel.hpp"
#include "./FChebM2LHandler.hpp"

class FTreeCoordinate;

/**
 * @author Matthias Messner(matthias.messner@inria.fr)
 * @class FChebKernel
 * @brief
 * Please read the license
 *
 * This kernels implement the Chebyshev 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
 */
43
template < class FReal, class CellClass, class ContainerClass,   class MatrixKernelClass, int ORDER, int NVALS = 1>
44
class FChebKernel
45
    : public FAbstractChebKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
46
{
47
    // private types
48
    typedef FChebM2LHandler<FReal, ORDER,MatrixKernelClass> M2LHandlerClass;
49

50
    // using from 
51
    typedef FAbstractChebKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
52
    AbstractBaseClass;
53

54 55
    /// Needed for P2P and M2L operators
    const MatrixKernelClass *const MatrixKernel;
56

57 58
    /// Needed for M2L operator
    FSmartPointer<  M2LHandlerClass,FSmartPointerMemory> M2LHandler;
59 60

public:
61 62 63 64 65
    /**
     * 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).
     */
66
    FChebKernel(const int inTreeHeight,  const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel,
67
                const FReal Epsilon)
68
    : FAbstractChebKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
69 70 71 72
                                                                                       inBoxWidth,
                                                                                       inBoxCenter),
      MatrixKernel(inMatrixKernel),
      M2LHandler(new M2LHandlerClass(MatrixKernel,Epsilon))
73 74 75 76 77
    {
        // read precomputed compressed m2l operators from binary file
        //M2LHandler->ReadFromBinaryFileAndSet();
        M2LHandler->ComputeAndCompressAndSet();
    }
78

79
    FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
80
        :   FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
81
    {
82

83
    }
84

85

86 87 88
    void P2M(CellClass* const LeafCell,
             const ContainerClass* const SourceParticles)
    {
89
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
90 91 92 93 94

        // 1) apply Sy
        AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                  LeafCell->getMultipole(0), SourceParticles);

95 96 97 98
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 2) apply B
            M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
        }
99
    }
100 101


102 103 104 105
    void M2M(CellClass* const FRestrict ParentCell,
             const CellClass*const FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
106 107 108 109 110 111 112 113 114 115 116
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 1) apply Sy
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
                if (ChildCells[ChildIndex]){
                    AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
                                                              ParentCell->getMultipole(idxRhs));
                }
            }
            // 2) apply B
            M2LHandler->applyB(ParentCell->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
        }
117 118 119 120
    }


//  void M2L(CellClass* const FRestrict TargetCell,
121
//                   const CellClass* SourceCells[],
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
//                   const int NumSourceCells,
//                   const int TreeLevel) const
//  {
//      const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
//      const FTreeCoordinate& cx = TargetCell->getCoordinate();
//      for (int idx=0; idx<NumSourceCells; ++idx) {
//          const FTreeCoordinate& cy = SourceCells[idx]->getCoordinate();
//          const int transfer[3] = {cy.getX()-cx.getX(),
//                                                           cy.getY()-cx.getY(),
//                                                           cy.getZ()-cx.getZ()};
//          M2LHandler->applyC(transfer, CellWidth,
//                                              SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
//                                              TargetCell->getLocal() + AbstractBaseClass::nnodes);
//      }
//  }

138 139
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int neighborPositions[], const int inSize, const int TreeLevel)  override {
140 141 142
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxRhs) + AbstractBaseClass::nnodes;
            const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
143 144 145 146
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idx = neighborPositions[idxExistingNeigh];
                M2LHandler->applyC(idx, CellWidth, SourceCells[idxExistingNeigh]->getMultipole(idxRhs) + AbstractBaseClass::nnodes,
                                   CompressedLocalExpansion);
147 148
            }
        }
149 150
    }

151 152
//  void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
//    const int neighborPositions[], const int inSize, const int TreeLevel)  override {
153 154 155
//      const unsigned int rank = M2LHandler.getRank();
//      FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
//      const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
156 157 158
//      for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
//          const int idx = neighborPositions[idxExistingNeigh];
//              FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idxExistingNeigh]->getMultipole())+AbstractBaseClass::nnodes,
159 160 161 162 163 164 165 166 167 168
//                                      MultipoleExpansion+idx*rank);
//      
//      M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
//  }


    void L2L(const CellClass* const FRestrict ParentCell,
             CellClass* FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
169 170 171
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 1) apply U
            M2LHandler->applyU(ParentCell->getLocal(idxRhs) + AbstractBaseClass::nnodes,
172
                               const_cast<CellClass*>(ParentCell)->getLocal(idxRhs));
173 174 175 176 177 178 179
            // 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));
                }
            }
        }
180
    }
181

182 183 184
    void L2P(const CellClass* const LeafCell,
             ContainerClass* const TargetParticles)
    {
185
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
186

187 188 189 190
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // 1) apply U
            M2LHandler->applyU(LeafCell->getLocal(idxRhs) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxRhs));
        }
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206

        //// 2.a) apply Sx
        //AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
        //                                          AbstractBaseClass::BoxWidthLeaf,
        //                                          LeafCell->getLocal(0),
        //                                          TargetParticles);
        //// 2.b) apply Px (grad Sx)
        //AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
        //                                                  AbstractBaseClass::BoxWidthLeaf,
        //                                                  LeafCell->getLocal(0),
        //                                                  TargetParticles);

        // 2.c) apply Sx and Px (grad Sx)
        AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                        LeafCell->getLocal(0), TargetParticles);

207
    }
208

209
    void P2P(const FTreeCoordinate& inPosition,
210
             ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict inSources,
211 212
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
213 214 215 216 217 218 219 220 221
        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);
        }
222 223 224 225 226 227
    }

    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
             ContainerClass* const FRestrict inTargets,
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
228 229 230 231 232 233
        int nbNeighborsToCompute = 0;
        while(nbNeighborsToCompute < inSize
              && neighborPositions[nbNeighborsToCompute] < 14){
            nbNeighborsToCompute += 1;
        }
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2P(inTargets,inNeighbors,nbNeighborsToCompute,MatrixKernel);
234 235 236 237
    }

    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
238
                   const ContainerClass* const inNeighbors[], const int /*neighborPositions*/[],
239 240
                   const int inSize) override {
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
241 242
    }

243 244 245 246 247 248
};


#endif //FCHEBKERNELS_HPP

// [--END--]