Commit b8583357 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Introduced the concept of extended bounding box in UniformTensorial kernel...

Introduced the concept of extended bounding box in UniformTensorial kernel (ONLY supported for NON_HOMOGENEOUS matrix kernels yet).
parent fef929b9
...@@ -59,8 +59,8 @@ protected: // PB for OptiDis ...@@ -59,8 +59,8 @@ protected: // PB for OptiDis
unsigned int node_ids[nnodes][3]; unsigned int node_ids[nnodes][3];
// 8 Non-leaf (i.e. M2M/L2L) interpolators // 8 Non-leaf (i.e. M2M/L2L) interpolators
// 1 per level if box is extended // x1 per level if box is extended
// TODO only 1 is required for all levels if extension is 0 // only 1 is required for all levels if extension is 0
FReal*** ChildParentInterpolator; FReal*** ChildParentInterpolator;
// Tree height (needed by M2M/L2L if cell width is extended) // Tree height (needed by M2M/L2L if cell width is extended)
......
...@@ -59,6 +59,9 @@ protected: ...@@ -59,6 +59,9 @@ protected:
const FReal BoxWidth; const FReal BoxWidth;
/// Width of a leaf cell box /// Width of a leaf cell box
const FReal BoxWidthLeaf; const FReal BoxWidthLeaf;
/// Extension of the box width ( same for all level! )
const FReal BoxWidthExtension;
/// Parameter to pass to matrix kernel (material specific or anything) /// Parameter to pass to matrix kernel (material specific or anything)
const double MatParam; const double MatParam;
...@@ -83,13 +86,17 @@ public: ...@@ -83,13 +86,17 @@ public:
FAbstractUnifKernel(const int inTreeHeight, FAbstractUnifKernel(const int inTreeHeight,
const FReal inBoxWidth, const FReal inBoxWidth,
const FPoint& inBoxCenter, const FPoint& inBoxCenter,
const FReal inBoxWidthExtension = 0.0,
const double inMatParam = 0.0) const double inMatParam = 0.0)
: Interpolator(new InterpolatorClass()), : Interpolator(new InterpolatorClass(inTreeHeight,
inBoxWidth,
inBoxWidthExtension)),
MatrixKernel(new MatrixKernelClass(inMatParam)), MatrixKernel(new MatrixKernelClass(inMatParam)),
TreeHeight(inTreeHeight), TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)), BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth), BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))), BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
BoxWidthExtension(inBoxWidthExtension),
MatParam(inMatParam) MatParam(inMatParam)
{ {
/* empty */ /* empty */
......
...@@ -49,7 +49,19 @@ class FUnifInterpolator : FNoCopyable ...@@ -49,7 +49,19 @@ class FUnifInterpolator : FNoCopyable
typedef FUnifTensor<ORDER> TensorType; typedef FUnifTensor<ORDER> TensorType;
unsigned int node_ids[nnodes][3]; unsigned int node_ids[nnodes][3];
FReal* ChildParentInterpolator[8];
// 8 Non-leaf (i.e. M2M/L2L) interpolators
// x1 per level if box is extended
// only 1 is required for all levels if extension is 0
FReal*** ChildParentInterpolator;
// Tree height (needed by M2M/L2L if cell width is extended)
const int TreeHeight;
// Root cell width (only used by M2M/L2L)
const FReal RootCellWidth;
// Cell width extension (only used by M2M/L2L, kernel handles extension for P2M/L2P)
const FReal CellWidthExtension;
// permutations (only needed in the tensor product interpolation case) // permutations (only needed in the tensor product interpolation case)
unsigned int perm[3][nnodes]; unsigned int perm[3][nnodes];
...@@ -95,29 +107,30 @@ class FUnifInterpolator : FNoCopyable ...@@ -95,29 +107,30 @@ class FUnifInterpolator : FNoCopyable
* S which is precomputed and reused for all M2M and L2L operations, ie for * S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations. * all non leaf inter/anterpolations.
*/ */
void initTensorM2MandL2L() void initTensorM2MandL2L(const int TreeLevel, const FReal ParentWidth)
{ {
FPoint ParentRoots[nnodes];
FReal ChildCoords[3][ORDER]; FReal ChildCoords[3][ORDER];
const FReal ParentWidth(2.);
const FPoint ParentCenter(0., 0., 0.);
FUnifTensor<ORDER>::setRoots(ParentCenter, ParentWidth, ParentRoots);
FPoint ChildCenter; FPoint ChildCenter;
const FReal ChildWidth(1.);
// Ratio of extended cell widths (definition: child ext / parent ext)
const FReal ExtendedCellRatio =
FReal(FReal(ParentWidth)/FReal(2.) + CellWidthExtension) / FReal(ParentWidth + CellWidthExtension);
// Child cell width
const FReal ChildWidth(2.*ExtendedCellRatio);
// loop: child cells // loop: child cells
for (unsigned int child=0; child<8; ++child) { for (unsigned int child=0; child<8; ++child) {
// set child info // set child info
FUnifTensor<ORDER>::setRelativeChildCenter(child, ChildCenter); FUnifTensor<ORDER>::setRelativeChildCenter(child, ChildCenter, ExtendedCellRatio);
FUnifTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords); FUnifTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory // allocate memory
ChildParentInterpolator[child] = new FReal [3 * ORDER*ORDER]; ChildParentInterpolator[TreeLevel][child] = new FReal [3 * ORDER*ORDER];
assembleInterpolator(ORDER, ChildCoords[0], ChildParentInterpolator[child]); assembleInterpolator(ORDER, ChildCoords[0], ChildParentInterpolator[TreeLevel][child]);
assembleInterpolator(ORDER, ChildCoords[1], ChildParentInterpolator[child] + 1 * ORDER*ORDER); assembleInterpolator(ORDER, ChildCoords[1], ChildParentInterpolator[TreeLevel][child] + 1 * ORDER*ORDER);
assembleInterpolator(ORDER, ChildCoords[2], ChildParentInterpolator[child] + 2 * ORDER*ORDER); assembleInterpolator(ORDER, ChildCoords[2], ChildParentInterpolator[TreeLevel][child] + 2 * ORDER*ORDER);
} }
...@@ -141,15 +154,51 @@ public: ...@@ -141,15 +154,51 @@ public:
/** /**
* Constructor: Initialize the Lagrange polynomials at the equispaced * Constructor: Initialize the Lagrange polynomials at the equispaced
* roots/interpolation point * roots/interpolation point
*
* PB: Input parameters ONLY affect the computation of the M2M/L2L ops.
* These parameters are ONLY required in the context of extended bbox.
* If no M2M/L2L is required then the interpolator can be built with
* the default ctor.
*/ */
explicit FUnifInterpolator() explicit FUnifInterpolator(const int inTreeHeight=3,
const FReal inRootCellWidth=FReal(1.),
const FReal inCellWidthExtension=FReal(0.))
: TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth),
CellWidthExtension(inCellWidthExtension)
{ {
// initialize root node ids // initialize root node ids
TensorType::setNodeIds(node_ids); TensorType::setNodeIds(node_ids);
// initialize interpolation operator for M2M and L2L (non leaf operations) // initialize interpolation operator for M2M and L2L (non leaf operations)
//this -> initM2MandL2L(); // non tensor-product interpolation
this -> initTensorM2MandL2L(); // tensor-product interpolation // allocate 8 arrays per level
ChildParentInterpolator = new FReal**[TreeHeight];
for (unsigned int l=0; l<TreeHeight; ++l){
ChildParentInterpolator[l] = new FReal*[8];
for (unsigned int c=0; c<8; ++c)
ChildParentInterpolator[l][c]=nullptr;
}
// Set number of non-leaf ios that actually need to be computed
unsigned int reducedTreeHeight; // = 2 + nb of computed nl ios
if(CellWidthExtension==0.) // if no cell extension, then ...
reducedTreeHeight = 3; // cmp only 1 non-leaf io
else
reducedTreeHeight = TreeHeight; // cmp 1 non-leaf io per level
// Init non-leaf interpolators
FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
CellWidth /= FReal(2.); // at level 2
for (unsigned int l=2; l<reducedTreeHeight; ++l) {
//this -> initM2MandL2L(l,CellWidth); // non tensor-product interpolation
this -> initTensorM2MandL2L(l,CellWidth); // tensor-product interpolation
// update cell width
CellWidth /= FReal(2.); // at level l+1
}
} }
...@@ -158,8 +207,10 @@ public: ...@@ -158,8 +207,10 @@ public:
*/ */
~FUnifInterpolator() ~FUnifInterpolator()
{ {
for (unsigned int child=0; child<8; ++child) for (unsigned int l=0; l<TreeHeight; ++l)
delete [] ChildParentInterpolator[child]; for (unsigned int child=0; child<8; ++child)
if(ChildParentInterpolator[l][child] != nullptr)
delete [] ChildParentInterpolator[l][child];
} }
...@@ -294,24 +345,25 @@ public: ...@@ -294,24 +345,25 @@ public:
// PB: Multidim version handled in kernel ! // PB: Multidim version handled in kernel !
void applyM2M(const unsigned int ChildIndex, void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion, const FReal *const ChildExpansion,
FReal *const ParentExpansion) const FReal *const ParentExpansion,
const unsigned int TreeLevel = 2) const
{ {
FReal Exp[nnodes], PermExp[nnodes]; FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER, ChildParentInterpolator[TreeLevel][ChildIndex], ORDER,
const_cast<FReal*>(ChildExpansion), ORDER, PermExp, ORDER); const_cast<FReal*>(ChildExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]]; for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 2 * ORDER*ORDER, ORDER, ChildParentInterpolator[TreeLevel][ChildIndex] + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER); Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]]; for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 1 * ORDER*ORDER, ORDER, ChildParentInterpolator[TreeLevel][ChildIndex] + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER); Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) ParentExpansion[perm[2][n]] += PermExp[n]; for (unsigned int n=0; n<nnodes; ++n) ParentExpansion[perm[2][n]] += PermExp[n];
...@@ -320,24 +372,25 @@ public: ...@@ -320,24 +372,25 @@ public:
void applyL2L(const unsigned int ChildIndex, void applyL2L(const unsigned int ChildIndex,
const FReal *const ParentExpansion, const FReal *const ParentExpansion,
FReal *const ChildExpansion) const FReal *const ChildExpansion,
const unsigned int TreeLevel = 2) const
{ {
FReal Exp[nnodes], PermExp[nnodes]; FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER, ChildParentInterpolator[TreeLevel][ChildIndex], ORDER,
const_cast<FReal*>(ParentExpansion), ORDER, PermExp, ORDER); const_cast<FReal*>(ParentExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]]; for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 2 * ORDER*ORDER, ORDER, ChildParentInterpolator[TreeLevel][ChildIndex] + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER); Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]]; for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-1) // ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.), FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 1 * ORDER*ORDER, ORDER, ChildParentInterpolator[TreeLevel][ChildIndex] + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER); Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) ChildExpansion[perm[2][n]] += PermExp[n]; for (unsigned int n=0; n<nnodes; ++n) ChildExpansion[perm[2][n]] += PermExp[n];
......
...@@ -87,11 +87,13 @@ public: ...@@ -87,11 +87,13 @@ public:
FUnifTensorialKernel(const int inTreeHeight, FUnifTensorialKernel(const int inTreeHeight,
const FReal inBoxWidth, const FReal inBoxWidth,
const FPoint& inBoxCenter, const FPoint& inBoxCenter,
const FReal inBoxWidthExtension,
const double inMatParam = 0.0) const double inMatParam = 0.0)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inMatParam), : FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(), M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight, inTreeHeight,
inBoxWidth) inBoxWidth,
inBoxWidthExtension)
{ } { }
...@@ -99,11 +101,13 @@ public: ...@@ -99,11 +101,13 @@ public:
const ContainerClass* const SourceParticles) const ContainerClass* const SourceParticles)
{ {
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension);
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
// 1) apply Sy // 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getMultipole(idxV*nRhs), SourceParticles); LeafCell->getMultipole(idxV*nRhs), SourceParticles);
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
...@@ -121,7 +125,7 @@ public: ...@@ -121,7 +125,7 @@ public:
void M2M(CellClass* const FRestrict ParentCell, void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells, const CellClass*const FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/) const int TreeLevel)
{ {
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
...@@ -132,8 +136,10 @@ public: ...@@ -132,8 +136,10 @@ public:
FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxMul)); FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxMul));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){ if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxMul), AbstractBaseClass::Interpolator->applyM2M(ChildIndex,
ParentCell->getMultipole(idxMul)); ChildCells[ChildIndex]->getMultipole(idxMul),
ParentCell->getMultipole(idxMul),
TreeLevel/*Cell width extension specific*/);
} }
} }
// 2) Apply Discete Fourier Transform // 2) Apply Discete Fourier Transform
...@@ -150,7 +156,8 @@ public: ...@@ -150,7 +156,8 @@ public:
const int TreeLevel) const int TreeLevel)
{ {
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel))); const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth)); const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(ExtendedCellWidth));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){ for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
...@@ -187,7 +194,7 @@ public: ...@@ -187,7 +194,7 @@ public:
void L2L(const CellClass* const FRestrict ParentCell, void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells, CellClass* FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/) const int TreeLevel)
{ {
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
...@@ -198,7 +205,10 @@ public: ...@@ -198,7 +205,10 @@ public:
// 2) apply Sx // 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){ if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxLoc), ChildCells[ChildIndex]->getLocal(idxLoc)); AbstractBaseClass::Interpolator->applyL2L(ChildIndex,
ParentCell->getLocal(idxLoc),
ChildCells[ChildIndex]->getLocal(idxLoc),
TreeLevel/*Cell width extension specific*/);
} }
} }
} }
...@@ -209,6 +219,8 @@ public: ...@@ -209,6 +219,8 @@ public:
ContainerClass* const TargetParticles) ContainerClass* const TargetParticles)
{ {
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension);
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
...@@ -220,11 +232,11 @@ public: ...@@ -220,11 +232,11 @@ public:
} }
// 2.a) apply Sx // 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getLocal(idxV*nLhs), TargetParticles); LeafCell->getLocal(idxV*nLhs), TargetParticles);
// 2.b) apply Px (grad Sx) // 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getLocal(idxV*nLhs), TargetParticles); LeafCell->getLocal(idxV*nLhs), TargetParticles);
}// NVALS }// NVALS
......
...@@ -38,6 +38,7 @@ PB: Compute() does not belong to the M2LHandler like it does in the Chebyshev ke ...@@ -38,6 +38,7 @@ PB: Compute() does not belong to the M2LHandler like it does in the Chebyshev ke
template <int ORDER, class MatrixKernelClass> template <int ORDER, class MatrixKernelClass>
static void Compute(const MatrixKernelClass *const MatrixKernel, static void Compute(const MatrixKernelClass *const MatrixKernel,
const FReal CellWidth, const FReal CellWidth,
const FReal CellWidthExtension,
FComplexe** &FC) FComplexe** &FC)
{ {
// dimensions of operators // dimensions of operators
...@@ -56,7 +57,8 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, ...@@ -56,7 +57,8 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
// interpolation points of source (Y) and target (X) cell // interpolation points of source (Y) and target (X) cell
FPoint X[nnodes], Y[nnodes]; FPoint X[nnodes], Y[nnodes];
// set roots of target cell (X) // set roots of target cell (X)
FUnifTensor<order>::setRoots(FPoint(0.,0.,0.), CellWidth, X); const FReal ExtendedCellWidth(CellWidth+CellWidthExtension);
FUnifTensor<order>::setRoots(FPoint(0.,0.,0.), ExtendedCellWidth, X);
// allocate memory and compute 316 m2l operators // allocate memory and compute 316 m2l operators
FReal** _C; FReal** _C;
...@@ -92,7 +94,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, ...@@ -92,7 +94,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
if (abs(i)>1 || abs(j)>1 || abs(k)>1) { if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
// set roots of source cell (Y) // set roots of source cell (Y)
const FPoint cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k)); const FPoint cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
FUnifTensor<order>::setRoots(cy, CellWidth, Y); FUnifTensor<order>::setRoots(cy, ExtendedCellWidth, Y);
// evaluate m2l operator // evaluate m2l operator
unsigned int ido=0; unsigned int ido=0;
for(unsigned int l=0; l<2*order-1; ++l) for(unsigned int l=0; l<2*order-1; ++l)
...@@ -193,6 +195,8 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl ...@@ -193,6 +195,8 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
/// M2L Operators (stored in Fourier space for each component of tensor) /// M2L Operators (stored in Fourier space for each component of tensor)
FSmartPointer< FComplexe*,FSmartArrayMemory> FC; FSmartPointer< FComplexe*,FSmartArrayMemory> FC;
const FReal CellWidthExtension; //<! extension of cells width
/// Utils /// Utils
typedef FUnifTensor<ORDER> TensorType; typedef FUnifTensor<ORDER> TensorType;
unsigned int node_diff[nnodes*nnodes]; unsigned int node_diff[nnodes*nnodes];
...@@ -215,8 +219,9 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl ...@@ -215,8 +219,9 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
public: public:
FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal) FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal inCellWidthExtension)
: opt_rc(rc/2+1), Dft() : CellWidthExtension(inCellWidthExtension),
Dft(), opt_rc(rc/2+1)
{ {
// init DFT // init DFT
const int steps[dimfft] = {rc}; const int steps[dimfft] = {rc};
...@@ -253,8 +258,11 @@ public: ...@@ -253,8 +258,11 @@ public:
if (FC[d]) throw std::runtime_error("M2L operator already set"); if (FC[d]) throw std::runtime_error("M2L operator already set");
// Compute matrix of interactions // Compute matrix of interactions
// reference cell width is arbitrarly set to 2.
// but it NEEDS to match the numerator of the scale factor in matrix kernel!
// Therefore box width extension is not yet supported for homog kernels
const FReal ReferenceCellWidth = FReal(2.); const FReal ReferenceCellWidth = FReal(2.);
Compute<order>(MatrixKernel,ReferenceCellWidth,FC); Compute<order>(MatrixKernel,ReferenceCellWidth, 0., FC);
// Compute memory usage // Compute memory usage
unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(FComplexe); unsigned long sizeM2L = 343*ncmp*opt_rc*sizeof(FComplexe);
...@@ -361,9 +369,11 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop ...@@ -361,9 +369,11 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
/// M2L Operators (stored in Fourier space for each component and each level) /// M2L Operators (stored in Fourier space for each component and each level)
FSmartPointer< FComplexe**,FSmartArrayMemory> FC; FSmartPointer< FComplexe**,FSmartArrayMemory> FC;
/// Homogeneity specific variables /// Homogeneity specific variables
const unsigned int TreeHeight; const unsigned int TreeHeight;
const FReal RootCellWidth; const FReal RootCellWidth;
const FReal CellWidthExtension; //<! extension of cells width
/// Utils /// Utils
typedef FUnifTensor<ORDER> TensorType; typedef FUnifTensor<ORDER> TensorType;
...@@ -387,10 +397,11 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop ...@@ -387,10 +397,11 @@ class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
public: public:
FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth) FUnifTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal inCellWidthExtension)
: TreeHeight(inTreeHeight), : TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth), RootCellWidth(inRootCellWidth),
opt_rc(rc/2+1), Dft() CellWidthExtension(inCellWidthExtension),
Dft(), opt_rc(rc/2+1)
{