FUnifKernel.hpp 9.11 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 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 47
// ===================================================================================
// 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".
// ===================================================================================
#ifndef FUNIFKERNEL_HPP
#define FUNIFKERNEL_HPP

#include "../../Utils/FGlobal.hpp"
#include "../../Utils/FTrace.hpp"
#include "../../Utils/FSmartPointer.hpp"

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

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
 */
template < class CellClass,	class ContainerClass,	class MatrixKernelClass, int ORDER, int NVALS = 1>
class FUnifKernel
  : public FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
  // private types
48
  typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
49 50 51 52 53

  // using from
  typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
  AbstractBaseClass;

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

57
  /// Needed for M2L operator
58 59
  const M2LHandlerClass M2LHandler;

60 61 62 63 64 65 66 67

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


  void P2M(CellClass* const LeafCell,
           const ContainerClass* const SourceParticles)
  {
    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
    for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
      // 1) apply Sy
      AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                LeafCell->getMultipole(idxRhs), SourceParticles);
      // 2) apply Discrete Fourier Transform
88
      M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxRhs), 
89 90 91 92 93 94 95 96 97 98 99 100
                                         LeafCell->getTransformedMultipole(idxRhs));

    }
  }


  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
101
      //FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
102 103 104 105 106 107 108
      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
109
      M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs), 
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
                                         ParentCell->getTransformedMultipole(idxRhs));
    }
  }


  //	void M2L(CellClass* const FRestrict TargetCell,
  //           const CellClass* SourceCells[343],
  //           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);
  //		}
  //	}

  void M2L(CellClass* const FRestrict TargetCell,
           const CellClass* SourceCells[343],
           const int /*NumSourceCells*/,
           const int TreeLevel)
  {
138
    const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
139
    const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
140

141 142 143 144 145
    for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
      FComplexe *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);

      for (int idx=0; idx<343; ++idx){
        if (SourceCells[idx]){
146
          M2LHandler.applyFC(idx, TreeLevel, scale, 
147
                              SourceCells[idx]->getTransformedMultipole(idxRhs),
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
                              TransformedLocalExpansion);
        }
      }
    }
  }

  //	void M2L(CellClass* const FRestrict TargetCell,
  //           const CellClass* SourceCells[343],
  //           const int NumSourceCells,
  //           const int TreeLevel) const
  //	{
  //		const unsigned int rank = M2LHandler.getRank();
  //		FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
  //		const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
  //		for (int idx=0; idx<343; ++idx)
  //			if (SourceCells[idx])
  //				FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idx]->getMultipole())+AbstractBaseClass::nnodes,
  //										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*/)
  {
    for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
      // 1) Apply Inverse Discete Fourier Transform
177
      M2LHandler.unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxRhs),
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
                                           const_cast<CellClass*>(ParentCell)->getLocal(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));
        }
      }
    }
  }

  void L2P(const CellClass* const LeafCell,
           ContainerClass* const TargetParticles)
  {
    const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));

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

      // 1)  Apply Inverse Discete Fourier Transform
196
      M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs), 
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
                                           const_cast<CellClass*>(LeafCell)->getLocal(idxRhs));

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

      // 2.b) apply Px (grad Sx)
      AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                        LeafCell->getLocal(idxRhs), 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 */)
  {
216
    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
217 218 219 220 221
  }


  void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                 ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
222 223
                 ContainerClass* const inNeighbors[27], const int /*inSize*/)
  {
224
    DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
225 226 227 228 229 230 231 232
  }

};


#endif //FUNIFKERNEL_HPP

// [--END--]