Commit f63e802b authored by Matthias Messner's avatar Matthias Messner

added ChebFMM for non-homogeneous kernel functions (eg. lenard-jones

potential); fixed some doxygen documentation
parent 3d293f9e
......@@ -117,7 +117,7 @@ public:
/**
*@brief Set the buffer capacity
*@param incapacity to change the capacity
*@param inCapacity to change the capacity
* If capacity given is lower than size elements after capacity are removed
*/
void setCapacity(const int inCapacity) {
......
......@@ -29,6 +29,7 @@ class FTreeCsvSaver {
public:
/** Constructor
* @param inBasefile is the output file name, you must put %d in it
* @param inIncludeHeader tells if header must be included
*/
FTreeCsvSaver(const char inBasefile[], const bool inIncludeHeader = false)
: includeHeader(inIncludeHeader), nbFrames(0) {
......
......@@ -27,9 +27,11 @@ class FTreeMpiCsvSaver {
static const int SizeOfLine = 56;
public:
/** Constructor
* @param inBasefile is the output file name, you must put %d in it
*/
/** Constructor
* @param inBasefile is the output file name, you must put %d in it
* @param communicator Mpi communicator
* @param inIncludeHeader tells if header must be included
*/
FTreeMpiCsvSaver(const char inBasefile[], const FMpi::FComm& communicator, const bool inIncludeHeader = false)
: comm(communicator.getComm()), includeHeader(inIncludeHeader), nbFrames(0) {
strcpy(basefile, inBasefile);
......
......@@ -11,13 +11,12 @@
#include "./FChebInterpolator.hpp"
class FTreeCoordinate;
template <KERNEL_FUNCCTION_IDENTIFIER Identifier> struct DirectInteactionComputer;
/**
* @author Matthias Messner(matthias.messner@inria.fr)
* @class FAbstractChebKernel
* @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.
......@@ -127,7 +126,7 @@ public:
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// only if target and source are not identical
this->directInteraction(Target, Source);
DirectInteactionComputer<MatrixKernelClass::Identifier>::compute(Target, Source);
// progress sources
iSources.gotoNext();
}
......@@ -140,7 +139,7 @@ public:
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// target and source cannot be identical
this->directInteraction(Target, Source);
DirectInteactionComputer<MatrixKernelClass::Identifier>::compute(Target, Source);
// progress sources
iSources.gotoNext();
}
......@@ -163,7 +162,7 @@ public:
while (iSources.hasNotFinished()) {
ParticleClass& Source = iSources.data();
// only if target and source are not identical
this->directInteractionMutual(Target, Source);
DirectInteactionComputer<MatrixKernelClass::Identifier>::computeMutual(Target, Source);
// progress sources
iSources.gotoNext();
}
......@@ -176,7 +175,7 @@ public:
while (iSources.hasNotFinished()) {
ParticleClass& Source = iSources.data();
// target and source cannot be identical
this->directInteractionMutual(Target, Source);
DirectInteactionComputer<MatrixKernelClass::Identifier>::computeMutual(Target, Source);
// progress sources
iSources.gotoNext();
}
......@@ -192,10 +191,15 @@ public:
}
};
private:
void directInteraction(ParticleClass& Target, const ParticleClass& Source) const // 34 overall flops
/*! Specialization for Laplace potential */
template <>
struct DirectInteactionComputer<ONE_OVER_R>
{
template <typename ParticleClass>
static void compute(ParticleClass& Target, const ParticleClass& Source) // 34 overall flops
{
FPoint xy(Source.getPosition() - Target.getPosition()); // 3 flops
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
......@@ -203,14 +207,17 @@ private:
xy.getZ()*xy.getZ()); // 1 + 15 + 5 = 21 flops
const FReal wt = Target.getPhysicalValue();
const FReal ws = Source.getPhysicalValue();
// potential
// laplace potential
Target.incPotential(one_over_r * ws); // 2 flops
// force
xy *= ((ws*wt) * (one_over_r*one_over_r*one_over_r)); // 5 flops
Target.incForces(xy.getX(), xy.getY(), xy.getZ()); // 3 flops
}
void directInteractionMutual(ParticleClass& Target, ParticleClass& Source) const // 39 overall flops
template <typename ParticleClass>
static void computeMutual(ParticleClass& Target, ParticleClass& Source) // 39 overall flops
{
FPoint xy(Source.getPosition() - Target.getPosition()); // 3 flops
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
......@@ -218,18 +225,70 @@ private:
xy.getZ()*xy.getZ()); // 1 + 15 + 5 = 21 flops
const FReal wt = Target.getPhysicalValue();
const FReal ws = Source.getPhysicalValue();
// potential
// laplace potential
Target.incPotential(one_over_r * ws); // 2 flops
Source.incPotential(one_over_r * wt); // 2 flops
// force
xy *= ((ws*wt) * (one_over_r*one_over_r*one_over_r)); // 5 flops
Target.incForces( xy.getX(), xy.getY(), xy.getZ()); // 3 flops
Source.incForces((-xy.getX()), (-xy.getY()), (-xy.getZ())); // 3 flops
}
};
/*! Specialization for Leonard-Jones potential */
template <>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL>
{
template <typename ParticleClass>
static void compute(ParticleClass& Target, const ParticleClass& Source) // 39 overall flops
{
FPoint xy(Source.getPosition() - Target.getPosition()); // 3 flops
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ()); // 1 + 15 + 5 = 21 flops
const FReal wt = Target.getPhysicalValue();
const FReal ws = Source.getPhysicalValue();
// lenard-jones potential
const FReal one_over_r3 = one_over_r * one_over_r * one_over_r;
const FReal one_over_r6 = one_over_r3 * one_over_r3;
Target.incPotential((one_over_r6*one_over_r6 - one_over_r6) * ws); // 2 flops
// force
const FReal one_over_r4 = one_over_r3 * one_over_r; // 1 flop
xy *= ((ws*wt) * (FReal(12.)*one_over_r6*one_over_r4*one_over_r4 - FReal(6.)*one_over_r4*one_over_r4)); // 9 flops
Target.incForces(xy.getX(), xy.getY(), xy.getZ()); // 3 flops
}
template <typename ParticleClass>
static void computeMutual(ParticleClass& Target, ParticleClass& Source) // 44 overall flops
{
FPoint xy(Source.getPosition() - Target.getPosition()); // 3 flops
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ()); // 1 + 15 + 5 = 21 flops
const FReal wt = Target.getPhysicalValue();
const FReal ws = Source.getPhysicalValue();
// lenard-jones potential
const FReal one_over_r3 = one_over_r * one_over_r * one_over_r;
const FReal one_over_r6 = one_over_r3 * one_over_r3;
Target.incPotential((one_over_r6*one_over_r6 - one_over_r6) * ws); // 2 flops
Source.incPotential((one_over_r6*one_over_r6 - one_over_r6) * wt); // 2 flops
// force
const FReal one_over_r4 = one_over_r3 * one_over_r; // 1 flop
xy *= ((ws*wt) * (FReal(12.)*one_over_r6*one_over_r4*one_over_r4 - FReal(6.)*one_over_r4*one_over_r4)); // 9 flops
Target.incForces( xy.getX(), xy.getY(), xy.getZ()); // 3 flops
Source.incForces((-xy.getX()), (-xy.getY()), (-xy.getZ())); // 3 flops
}
};
#endif //FCHEBKERNELS_HPP
// [--END--]
......@@ -6,9 +6,13 @@
#include "../../Utils/FNoCopyable.hpp"
#include "../../Utils/FMath.hpp"
// extendable
enum KERNEL_FUNCCTION_IDENTIFIER {ONE_OVER_R,
ONE_OVER_R_SQUARED,
LEONARD_JONES_POTENTIAL};
enum {ONE_OVER_R = 1,
ONE_OVER_R_SQUARED = 2};
// probably not extedable :)
enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
/**
......@@ -19,14 +23,21 @@ enum {ONE_OVER_R = 1,
struct FChebAbstractMatrixKernel : FNoCopyable
{
virtual FReal evaluate(const FPoint&, const FPoint&) const = 0;
virtual FReal getScaleFactor(FReal) const = 0;
// I need both functions because required arguments are not always given
virtual FReal getScaleFactor(const FReal, const int) const = 0;
virtual FReal getScaleFactor(const FReal) const = 0;
};
/// One over r
struct FChebMatrixKernelR : FChebAbstractMatrixKernel
{
enum {Identifier = ONE_OVER_R};
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R;
FChebMatrixKernelR() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -35,15 +46,28 @@ struct FChebMatrixKernelR : FChebAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
FReal getScaleFactor(FReal CellWidth) const
{ return FReal(2.) / CellWidth; }
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(2.) / CellWidth;
}
};
/// One over r^2
struct FChebMatrixKernelRR : FChebAbstractMatrixKernel
{
enum {Identifier = ONE_OVER_R_SQUARED};
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
FChebMatrixKernelRR() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -52,11 +76,56 @@ struct FChebMatrixKernelRR : FChebAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
FReal getScaleFactor(FReal CellWidth) const
{ return FReal(4.) / (CellWidth*CellWidth); }
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(4.) / CellWidth;
}
};
/// One over r^12 - One over r^6
struct FChebMatrixKernelLJ : FChebAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = LEONARD_JONES_POTENTIAL;
FChebMatrixKernelLJ() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
const FReal r = xy.norm();
const FReal r3 = r*r*r;
const FReal one_over_r6 = FReal(1.) / (r3*r3);
//return one_over_r6 * one_over_r6;
//return one_over_r6;
return one_over_r6 * one_over_r6 - one_over_r6;
}
FReal getScaleFactor(const FReal, const int) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
return FReal(1.);
}
FReal getScaleFactor(const FReal) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
return FReal(1.);
}
};
#endif // FCHEBMATRIXKERNEL_HPP
// [--END--]
......@@ -36,7 +36,7 @@ class FChebSymKernel
: public FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER>
{
typedef FAbstractChebKernel<ParticleClass, CellClass, ContainerClass, MatrixKernelClass, ORDER> AbstractBaseClass;
typedef SymmetryHandler<ORDER> SymmetryHandlerClass;
typedef SymmetryHandler<ORDER, MatrixKernelClass::Type> SymmetryHandlerClass;
enum {nnodes = AbstractBaseClass::nnodes};
/// Needed for handling all symmetries
......@@ -88,7 +88,7 @@ public:
const FReal inBoxWidth,
const FReal Epsilon)
: AbstractBaseClass(inTreeHeight, inBoxCenter, inBoxWidth),
SymHandler(new SymmetryHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon)),
SymHandler(new SymmetryHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon, inBoxWidth, inTreeHeight)),
Loc(NULL), Mul(NULL), countExp(NULL)
{
this->allocateMemoryForPermutedExpansions();
......@@ -153,7 +153,6 @@ public:
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
......@@ -175,18 +174,19 @@ public:
// multiply (mat-mat-mul)
FReal Compressed [nnodes * 24];
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel->getScaleFactor(CellWidth));
const FReal scale = AbstractBaseClass::MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
for (unsigned int pidx=0; pidx<343; ++pidx) {
const unsigned int count = countExp[pidx];
if (count) {
const unsigned int rank = SymHandler->LowRank[pidx];
const unsigned int rank = SymHandler->getLowRank(TreeLevel, pidx);
// rank * count * (2*nnodes-1) flops
FBlas::gemtm(nnodes, rank, count, FReal(1.),
SymHandler->K[pidx]+rank*nnodes, nnodes, Mul[pidx], nnodes, Compressed, rank);
const_cast<FReal*>(SymHandler->getK(TreeLevel, pidx))+rank*nnodes,
nnodes, Mul[pidx], nnodes, Compressed, rank);
// nnodes *count * (2*rank-1) flops
FBlas::gemm( nnodes, rank, count, scale,
SymHandler->K[pidx], nnodes, Compressed, rank, Loc[pidx], nnodes);
const_cast<FReal*>(SymHandler->getK(TreeLevel, pidx)),
nnodes, Compressed, rank, Loc[pidx], nnodes);
}
}
......@@ -207,6 +207,84 @@ public:
}
// template <KERNEL_FUNCTION_TYPE TYPE>
// struct BlockedOuterProduct;
//
// template<>
// struct BlockedOuterProduct<HOMOGENEOUS>
// {
// static void apply() const
// {
// FReal Compressed [nnodes * 24];
// 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) {
// const unsigned int rank = SymHandler->LowRank[pidx];
// // rank * count * (2*nnodes-1) flops
// FBlas::gemtm(nnodes, rank, count, FReal(1.),
// SymHandler->K[pidx]+rank*nnodes, nnodes, Mul[pidx], nnodes, Compressed, rank);
// // nnodes *count * (2*rank-1) flops
// FBlas::gemm( nnodes, rank, count, scale,
// SymHandler->K[pidx], nnodes, Compressed, rank, Loc[pidx], nnodes);
// }
// }
// }
// };
//template <class ParticleClass, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
//template<>
//void FChebSymKernel::BlockedOuterProduct<NON_HOMOGENEOUS>()
//{
// FReal Compressed [nnodes * 24];
// for (unsigned int pidx=0; pidx<343; ++pidx) {
// const unsigned int count = countExp[pidx];
// if (count) {
// const unsigned int rank = SymHandler->LowRank[pidx];
// // rank * count * (2*nnodes-1) flops
// FBlas::gemtm(nnodes, rank, count, FReal(1.),
// SymHandler->K[pidx]+rank*nnodes, nnodes, Mul[pidx], nnodes, Compressed, rank);
// // nnodes *count * (2*rank-1) flops
// FBlas::gemm( nnodes, rank, count, FReal(1.),
// SymHandler->K[pidx], nnodes, Compressed, rank, Loc[pidx], nnodes);
// }
// }
//}
//template <KERNEL_FUNCTION_TYPE TYPE>
//struct blocked_outer_product { static void apply(const int TreeLevel); };
//
//template <>
//struct blocked_outer_product<HOMOGENEOUS>
//{
// static void apply(const int TreeLevel)
// {
// FReal Compressed [nnodes * 24];
// 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) {
// const unsigned int rank = SymHandler->LowRank[pidx];
// // rank * count * (2*nnodes-1) flops
// FBlas::gemtm(nnodes, rank, count, FReal(1.),
// SymHandler->K[pidx]+rank*nnodes, nnodes, Mul[pidx], nnodes, Compressed, rank);
// // nnodes *count * (2*rank-1) flops
// FBlas::gemm( nnodes, rank, count, scale,
// SymHandler->K[pidx], nnodes, Compressed, rank, Loc[pidx], nnodes);
// }
// }
// }
//};
//
//template <>
//struct blocked_outer_product<NON_HOMOGENEOUS>
//{
// static void apply(const int TreeLevel) {}
//};
/*
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
......@@ -345,19 +423,6 @@ public:
#endif //FCHEBSYMKERNELS_HPP
......
......@@ -119,17 +119,17 @@ static void fACA(FReal *const K, const double eps, FReal* &U, FReal* &V, unsigne
FACASVD is defined or not, either ACA+SVD or only SVD is used to compress
them. */
template <int ORDER, typename MatrixKernelClass>
static void precompute(const MatrixKernelClass *const MatrixKernel, const double Epsilon, FReal* K[343], int LowRank[343])
static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal CellWidth,
const double Epsilon, FReal* K[343], int LowRank[343])
{
static const unsigned int nnodes = ORDER*ORDER*ORDER;
std::cout << "\nComputing 16 far-field interactions for cells of width w = " << CellWidth << std::endl;
// set all to zero
for (unsigned int t=0; t<343; ++t) { K[t] = NULL; LowRank[t] = 0; }
static const unsigned int nnodes = ORDER*ORDER*ORDER;
// interpolation points of source (Y) and target (X) cell
FPoint X[nnodes], Y[nnodes];
// set roots of target cell (X)
FChebTensor<ORDER>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
FChebTensor<ORDER>::setRoots(FPoint(0.,0.,0.), CellWidth, X);
// temporary matrix
FReal* U = new FReal [nnodes*nnodes];
......@@ -146,8 +146,8 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const double
for (int k=0; k<=j; ++k) {
// assemble matrix
const FPoint cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<ORDER>::setRoots(cy, FReal(2.), Y);
const FPoint cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
FChebTensor<ORDER>::setRoots(cy, CellWidth, Y);
for (unsigned int n=0; n<nnodes; ++n)
for (unsigned int m=0; m<nnodes; ++m)
U[n*nnodes + m] = MatrixKernel->evaluate(X[m], Y[n]);
......@@ -206,7 +206,6 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const double
INFO = FBlas::gesvd(aca_rank, aca_rank, phi, S, VT, aca_rank, LWORK, WORK);
if (INFO!=0) throw std::runtime_error("SVD did not converge with " + INFO);
rank = getRank(S, aca_rank, Epsilon);
//std::cout << " - recompression with SVD leads to rank = " << rank << std::endl;
}
const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
......@@ -310,15 +309,19 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const double
because it allows us to use to associate the far-field interactions based on
the index \f$t = 7^2(i+3) + 7(j+3) + (k+3)\f$ where \f$(i,j,k)\f$ denotes
the relative position of the source cell to the target cell. */
template <int ORDER, KERNEL_FUNCTION_TYPE TYPE> class SymmetryHandler;
/*! Specialization for homogeneous kernel functions */
template <int ORDER>
class SymmetryHandler
class SymmetryHandler<ORDER, HOMOGENEOUS>
{
static const unsigned int nnodes = ORDER*ORDER*ORDER;
public:
// M2L operators
FReal* K[343];
int LowRank[343];
public:
// permutation vectors and permutated indices
unsigned int pvectors[343][nnodes];
......@@ -327,7 +330,8 @@ public:
/** Constructor: with 16 small SVDs */
template <typename MatrixKernelClass>
SymmetryHandler(const MatrixKernelClass *const MatrixKernel, const double Epsilon)
SymmetryHandler(const MatrixKernelClass *const MatrixKernel, const double Epsilon,
const FReal, const unsigned int)
{
// init all 343 item to zero, because effectively only 16 exist
for (unsigned int t=0; t<343; ++t) {
......@@ -347,7 +351,8 @@ public:
}
// precompute 16 M2L operators
precompute<ORDER>(MatrixKernel, Epsilon, K, LowRank);
const FReal ReferenceCellWidth = FReal(2.);
precompute<ORDER>(MatrixKernel, ReferenceCellWidth, Epsilon, K, LowRank);
}
......@@ -358,6 +363,107 @@ public:
for (unsigned int t=0; t<343; ++t) if (K[t]!=NULL) delete [] K[t];
}
/*! return the t-th approximated far-field interactions*/
const FReal *const getK(const unsigned int, const unsigned int t) const
{ return K[t]; }
/*! return the t-th approximated far-field interactions*/
const int getLowRank(const unsigned int, const unsigned int t) const
{ return LowRank[t]; }
};
/*! Specialization for non-homogeneous kernel functions */
template <int ORDER>
class SymmetryHandler<ORDER, NON_HOMOGENEOUS>
{
static const unsigned int nnodes = ORDER*ORDER*ORDER;
// Height of octree; needed only in the case of non-homogeneous kernel functions
const unsigned int TreeHeight;
// M2L operators for all levels in the octree
FReal*** K;
int** LowRank;
public:
// permutation vectors and permutated indices
unsigned int pvectors[343][nnodes];
unsigned int pindices[343];
/** Constructor: with 16 small SVDs */
template <typename MatrixKernelClass>
SymmetryHandler(const MatrixKernelClass *const MatrixKernel, const double Epsilon,
const FReal RootCellWidth, const unsigned int inTreeHeight)
: TreeHeight(inTreeHeight)
{
// init all 343 item to zero, because effectively only 16 exist
K = new FReal** [TreeHeight];
LowRank = new int* [TreeHeight];
K[0] = NULL; K[1] = NULL;
LowRank[0] = NULL; LowRank[1] = NULL;
for (unsigned int l=2; l<TreeHeight; ++l) {
K[l] = new FReal* [343];
LowRank[l] = new int [343];
for (unsigned int t=0; t<343; ++t) {
K[l][t] = NULL;
LowRank[l][t] = 0;
}
}
// set permutation vector and indices
const FChebSymmetries<ORDER> Symmetries;
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
const unsigned int idx = ((i+3) * 7 + (j+3)) * 7 + (k+3);
pindices[idx] = 0;
if (abs(i)>1 || abs(j)>1 || abs(k)>1)
pindices[idx] = Symmetries.getPermutationArrayAndIndex(i,j,k, pvectors[idx]);
}
// precompute 16 M2L operators at all levels having far-field interactions
FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
CellWidth /= FReal(2.); // at level 2
for (unsigned int l=2; l<TreeHeight; ++l) {
precompute<ORDER>(MatrixKernel, CellWidth, Epsilon, K[l], LowRank[l]);
CellWidth /= FReal(2.); // at level l+1
}
}
/** Destructor */
~SymmetryHandler()
{
for (unsigned int l=0; l<TreeHeight; ++l) {
if (K[l]!=NULL) {
for (unsigned int t=0; t<343; ++t) if (K[l][t]!=NULL) delete [] K[l][t];
delete [] K[l];
}
if (LowRank[l]!=NULL) delete [] LowRank[l];
}
delete [] K;
delete [] LowRank;
}
/*! return the t-th approximated far-field interactions*/
const FReal *const getK(const unsigned int l, const unsigned int t) const
{ return K[l][t]; }
/*! return the t-th approximated far-field interactions*/
const int getLowRank(const unsigned int l, const unsigned int t) const
{ return LowRank[l][t]; }
};