Commit 8b54072f authored by Matthias Messner's avatar Matthias Messner
Browse files

added abstract base class FAbstractChebKernel, changed name of ChebKernels to ChebKernel

parent ef6fd9f8
#ifndef FCHEBKERNELS_HPP
#define FCHEBKERNELS_HPP
#ifndef FABSTRACTCHEBKERNEL_HPP
#define FABSTRACTCHEBKERNEL_HPP
// [--License--]
#include "../Utils/FGlobal.hpp"
......@@ -9,13 +9,12 @@
#include "../Components/FAbstractKernels.hpp"
#include "./FChebInterpolator.hpp"
#include "./FChebM2LHandler.hpp"
class FTreeCoordinate;
/**
* @author Matthias Messner(matthias.messner@inria.fr)
* @class FChebKernels
* @class FAbstractChebKernel
* @brief
* Please read the license
*
......@@ -30,16 +29,14 @@ class FTreeCoordinate;
* @tparam ORDER Chebyshev interpolation order
*/
template <class ParticleClass, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
class FChebKernels : public FAbstractKernels<ParticleClass, CellClass, ContainerClass>
class FAbstractChebKernel : public FAbstractKernels<ParticleClass, CellClass, ContainerClass>
{
protected:
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebInterpolator<ORDER> InterpolatorClass;
typedef FChebM2LHandler<ORDER,MatrixKernelClass> M2LHandlerClass;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for M2L operator
FSmartPointer< M2LHandlerClass,FSmartPointerMemory> M2LHandler;
/// Needed for P2P operator
const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Height of the entire oct-tree
......@@ -50,8 +47,6 @@ class FChebKernels : public FAbstractKernels<ParticleClass, CellClass, Container
const FReal BoxWidth;
/// Width of a leaf cell box
const FReal BoxWidthLeaf;
/// Prescribed accuracy of the compressed M2L operators
const FReal Epsilon;
/**
* Compute center of leaf cell from its tree coordinate.
......@@ -65,147 +60,51 @@ class FChebKernels : public FAbstractKernels<ParticleClass, CellClass, Container
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
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).
*/
FChebKernels(const int inTreeHeight,
const F3DPosition& inBoxCenter,
const FReal inBoxWidth,
const FReal inEpsilon)
FAbstractChebKernel(const int inTreeHeight,
const F3DPosition& inBoxCenter,
const FReal inBoxWidth)
: Interpolator(new InterpolatorClass()),
M2LHandler(new M2LHandlerClass(inEpsilon)),
MatrixKernel(new MatrixKernelClass()),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
Epsilon(inEpsilon)
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
{
// read precomputed compressed m2l operators from binary file
M2LHandler->ReadFromBinaryFileAndSet();
//M2LHandler->ComputeAndCompressAndSet();
/* empty */
}
void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles)
{
// 1) apply Sy
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
Interpolator->applyP2M(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getMultipole(),
SourceParticles);
// 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(),
LeafCell->getMultipole() + nnodes);
}
virtual void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles) = 0;
void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel)
{
// 1) apply Sy
FBlas::scal(nnodes*2, FReal(0.), ParentCell->getMultipole());
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(),
ParentCell->getMultipole());
// 2) apply B
M2LHandler->applyB(ParentCell->getMultipole(),
ParentCell->getMultipole() + nnodes);
}
virtual void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
// 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() + nnodes,
// TargetCell->getLocal() + nnodes);
// }
// }
virtual void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int TreeLevel) = 0;
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int TreeLevel)
{
FReal *const CompressedLocalExpansion = TargetCell->getLocal() + nnodes;
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for (int idx=0; idx<343; ++idx)
if (SourceCells[idx])
M2LHandler->applyC(idx, CellWidth,
SourceCells[idx]->getMultipole() + nnodes,
CompressedLocalExpansion);
}
// 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())+nnodes,
// MultipoleExpansion+idx*rank);
//
// M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + nnodes);
// }
virtual void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int TreeLevel)
{
// 1) apply U
M2LHandler->applyU(ParentCell->getLocal() + nnodes,
const_cast<CellClass*>(ParentCell)->getLocal());
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
Interpolator->applyL2L(ChildIndex,
ParentCell->getLocal(),
ChildCells[ChildIndex]->getLocal());
}
void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles)
{
// 1) apply U
M2LHandler->applyU(LeafCell->getLocal() + nnodes,
const_cast<CellClass*>(LeafCell)->getLocal());
// 2.a) apply Sx
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
Interpolator->applyL2P(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// 2.b) apply Px (grad Sx)
Interpolator->applyL2PGradient(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
}
virtual void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles) = 0;
void P2P(const FTreeCoordinate& LeafCellCoordinate,
ContainerClass* const FRestrict TargetParticles,
......
#ifndef FCHEBKERNEL_HPP
#define FCHEBKERNEL_HPP
// [--License--]
#include "../Utils/FGlobal.hpp"
#include "../Utils/FTrace.hpp"
#include "../Utils/FSmartPointer.hpp"
#include "./FAbstractChebKernel.hpp"
//#include "./FChebInterpolator.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 ParticleClass Type of particle
* @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 ParticleClass, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
class FChebKernel
: public FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>
{
// private types
typedef FChebM2LHandler<ORDER,MatrixKernelClass> M2LHandlerClass;
// using from
typedef FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>
AbstractBaseClass;
/// Needed for M2L operator
FSmartPointer< M2LHandlerClass,FSmartPointerMemory> M2LHandler;
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).
*/
FChebKernel(const int inTreeHeight,
const F3DPosition& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>(inTreeHeight,
inBoxCenter,
inBoxWidth),
M2LHandler(new M2LHandlerClass(Epsilon))
{
// read precomputed compressed m2l operators from binary file
M2LHandler->ReadFromBinaryFileAndSet();
//M2LHandler->ComputeAndCompressAndSet();
}
void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles)
{
// 1) apply Sy
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(),
SourceParticles);
// 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(),
LeafCell->getMultipole() + AbstractBaseClass::nnodes);
}
void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel)
{
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole());
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
AbstractBaseClass::Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(),
ParentCell->getMultipole());
// 2) apply B
M2LHandler->applyB(ParentCell->getMultipole(),
ParentCell->getMultipole() + AbstractBaseClass::nnodes);
}
// 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)
{
FReal *const CompressedLocalExpansion = TargetCell->getLocal() + AbstractBaseClass::nnodes;
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for (int idx=0; idx<343; ++idx)
if (SourceCells[idx])
M2LHandler->applyC(idx, CellWidth,
SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
CompressedLocalExpansion);
}
// 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)
{
// 1) apply U
M2LHandler->applyU(ParentCell->getLocal() + AbstractBaseClass::nnodes,
const_cast<CellClass*>(ParentCell)->getLocal());
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
AbstractBaseClass::Interpolator->applyL2L(ChildIndex,
ParentCell->getLocal(),
ChildCells[ChildIndex]->getLocal());
}
void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles)
{
// 1) apply U
M2LHandler->applyU(LeafCell->getLocal() + AbstractBaseClass::nnodes,
const_cast<CellClass*>(LeafCell)->getLocal());
// 2.a) apply Sx
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
}
};
#endif //FCHEBKERNELS_HPP
// [--END--]
#ifndef FCHEBSYMKERNELS_HPP
#define FCHEBSYMKERNELS_HPP
#ifndef FCHEBSYMKERNEL_HPP
#define FCHEBSYMKERNEL_HPP
// [--License--]
#include "../Utils/FGlobal.hpp"
#include "../Utils/FTrace.hpp"
#include "../Utils/FSmartPointer.hpp"
#include "../Components/FAbstractKernels.hpp"
#include "./FAbstractChebKernel.hpp"
#include "./FChebInterpolator.hpp"
#include "./FChebSymmetries.hpp"
......@@ -17,7 +17,7 @@ class FTreeCoordinate;
/**
* @author Matthias Messner(matthias.messner@inria.fr)
* @class FChebSymKernels
* @class FChebSymKernel
* @brief
* Please read the license
*
......@@ -33,49 +33,24 @@ class FTreeCoordinate;
* @tparam ORDER Chebyshev interpolation order
*/
template <class ParticleClass, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
class FChebSymKernels : public FAbstractKernels<ParticleClass, CellClass, ContainerClass>
class FChebSymKernel
: public FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>
{
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebInterpolator<ORDER> InterpolatorClass;
typedef FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>
AbstractBaseClass;
enum {nnodes = AbstractBaseClass::nnodes};
/// Handler to deal with all symmetries: Stores permutation indices and
/// vectors to reduce 343 different interactions to 16 only.
struct SymmetryHandler;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for P2P operator
const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Needed for handling all symmetries
const FSmartPointer<SymmetryHandler,FSmartPointerMemory> SymHandler;
/// Height of the entire oct-tree
const unsigned int TreeHeight;
/// Corner of oct-tree box
const F3DPosition BoxCorner;
/// Width of oct-tree box
const FReal BoxWidth;
/// Width of a leaf cell box
const FReal BoxWidthLeaf;
// permuted local and multipole expansions
FReal* Loc[343];
FReal* Mul[343];
unsigned int countExp[343];
/**
* Compute center of leaf cell from its tree coordinate.
* @param[in] Coordinate tree coordinate
* @return center of leaf cell
*/
const F3DPosition getLeafCellCenter(const FTreeCoordinate& Coordinate) const
{
return F3DPosition(BoxCorner.getX() + (FReal(Coordinate.getX()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getY() + (FReal(Coordinate.getY()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
FReal** Loc;
FReal** Mul;
unsigned int* countExp;
......@@ -84,11 +59,22 @@ class FChebSymKernels : public FAbstractKernels<ParticleClass, CellClass, Contai
*/
void allocateMemoryForPermutedExpansions()
{
assert(Loc==NULL && Mul==NULL && countExp==NULL);
Loc = new FReal* [343];
Mul = new FReal* [343];
countExp = new unsigned int [343];
// set all 343 to NULL
for (unsigned int idx=0; idx<343; ++idx) {
Mul[idx] = Loc[idx] = NULL;
}
// init only 16 of 343 possible translations due to symmetries
for (int i=2; i<=3; ++i)
for (int j=0; j<=i; ++j)
for (int k=0; k<=j; ++k) {
const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
assert(Mul[idx]!=NULL || Loc[idx]!=NULL);
assert(Mul[idx]==NULL || Loc[idx]==NULL);
Mul[idx] = new FReal [30 * nnodes];
Loc[idx] = new FReal [30 * nnodes];
}
......@@ -102,32 +88,24 @@ public:
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebSymKernels(const int inTreeHeight,
const F3DPosition& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: Interpolator(new InterpolatorClass()),
MatrixKernel(new MatrixKernelClass()),
SymHandler(new SymmetryHandler(MatrixKernel.getPtr(), Epsilon)),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
FChebSymKernel(const int inTreeHeight,
const F3DPosition& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: AbstractBaseClass(inTreeHeight, inBoxCenter, inBoxWidth),
SymHandler(new SymmetryHandler(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon)),
Loc(NULL), Mul(NULL), countExp(NULL)
{
this->allocateMemoryForPermutedExpansions();
}
/** Copy constructor */
FChebSymKernels(const FChebSymKernels& other)
: Interpolator(other.Interpolator),
MatrixKernel(other.MatrixKernel),
FChebSymKernel(const FChebSymKernel& other)
: AbstractBaseClass(other),
SymHandler(other.SymHandler),
TreeHeight(other.TreeHeight),
BoxCorner(other.BoxCorner),
BoxWidth(other.BoxWidth),
BoxWidthLeaf(other.BoxWidthLeaf)
Loc(NULL), Mul(NULL), countExp(NULL)
{
this->allocateMemoryForPermutedExpansions();
}
......@@ -135,12 +113,15 @@ public:
/** Destructor */
~FChebSymKernels()
~FChebSymKernel()
{
for (unsigned int t=0; t<343; ++t) {
if (Loc[t]!=NULL) delete [] Loc[t];
if (Mul[t]!=NULL) delete [] Mul[t];
}
if (Loc!=NULL) delete [] Loc;
if (Mul!=NULL) delete [] Mul;
if (countExp!=NULL) delete [] countExp;
}
......@@ -150,10 +131,10 @@ public:
{
// apply Sy
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
Interpolator->applyP2M(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getMultipole(),
SourceParticles);
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(),
SourceParticles);
}
......@@ -166,9 +147,9 @@ public:
FBlas::scal(nnodes*2, FReal(0.), ParentCell->getMultipole());
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(),
ParentCell->getMultipole());
AbstractBaseClass::Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(),
ParentCell->getMultipole());
}
......@@ -193,8 +174,8 @@ public:
// multiply (mat-mat-mul)
FReal Compressed [nnodes * 30];
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel->getScaleFactor(CellWidth));
for (unsigned int pidx=0; pidx<343; ++pidx) {
const unsigned int count = countExp[pidx];
if (count) {
......@@ -222,11 +203,11 @@ public:
}
/*
void M2L(CellClass* const FRestrict TargetCell,