Commit 370adb44 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

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

Introduced the concept of extended bounding box in ChebyshevTensorial kernel (ONLY supported for NON_HOMOGENEOUS kernels yet) and provided corresponding tests for interpolator.
parent 45834df4
......@@ -59,6 +59,9 @@ protected:
const FReal BoxWidth;
/// Width of a leaf cell box
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)
const double MatParam;
......@@ -83,13 +86,17 @@ public:
FAbstractChebKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal inBoxWidthExtension = 0.0,
const double inMatParam = 0.0)
: Interpolator(new InterpolatorClass()),
: Interpolator(new InterpolatorClass(inTreeHeight,
inBoxWidth,
inBoxWidthExtension)),
MatrixKernel(new MatrixKernelClass(inMatParam)),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
BoxWidthExtension(inBoxWidthExtension),
MatParam(inMatParam)
{
/* empty */
......
......@@ -56,8 +56,20 @@ protected: // PB for OptiDis
FReal T_of_roots[ORDER][ORDER];
FReal T[ORDER * (ORDER-1)];
unsigned int node_ids[nnodes][3];
FReal* ChildParentInterpolator[8];
unsigned int node_ids[nnodes][3];
// 8 Non-leaf (i.e. M2M/L2L) interpolators
// 1 per level if box is extended
// TODO 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)
unsigned int perm[3][nnodes];
......@@ -259,85 +271,103 @@ protected: // PB for OptiDis
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations.
* This is a sub-optimal version : complexity p^6.
* This function handles extended cells.
*/
void initM2MandL2L()
{
FPoint ParentRoots[nnodes], ChildRoots[nnodes];
const FReal ParentWidth(2.);
const FPoint ParentCenter(0., 0., 0.);
FChebTensor<ORDER>::setRoots(ParentCenter, ParentWidth, ParentRoots);
void initM2MandL2L(const int TreeLevel, const FReal ParentWidth)
{
FPoint ChildRoots[nnodes];
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);
// loop: child cells
for (unsigned int child=0; child<8; ++child) {
// Child cell center and width
FPoint ChildCenter;
const FReal ChildWidth(2.*ExtendedCellRatio);
// allocate memory
ChildParentInterpolator[child] = new FReal [nnodes * nnodes];
// loop: child cells
for (unsigned int child=0; child<8; ++child) {
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
FChebTensor<ORDER>::setRoots(ChildCenter, ChildWidth, ChildRoots);
// allocate memory
ChildParentInterpolator[TreeLevel][child] = new FReal [nnodes * nnodes];
// assemble child - parent - interpolator
assembleInterpolator(nnodes, ChildRoots, ChildParentInterpolator[child]);
}
}
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter, ExtendedCellRatio);
FChebTensor<ORDER>::setRoots(ChildCenter, ChildWidth, ChildRoots);
/**
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations.
*/
void initTensorM2MandL2L()
{
FPoint ParentRoots[nnodes];
FReal ChildCoords[3][ORDER];
const FReal ParentWidth(2.);
const FPoint ParentCenter(0., 0., 0.);
FChebTensor<ORDER>::setRoots(ParentCenter, ParentWidth, ParentRoots);
FPoint ChildCenter;
const FReal ChildWidth(1.);
// loop: child cells
for (unsigned int child=0; child<8; ++child) {
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
FChebTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory
ChildParentInterpolator[child] = new FReal [3 * ORDER*ORDER];
assembleInterpolator(ORDER, ChildCoords[0], ChildParentInterpolator[child]);
assembleInterpolator(ORDER, ChildCoords[1], ChildParentInterpolator[child] + 1 * ORDER*ORDER);
assembleInterpolator(ORDER, ChildCoords[2], ChildParentInterpolator[child] + 2 * ORDER*ORDER);
}
// assemble child - parent - interpolator
assembleInterpolator(nnodes, ChildRoots, ChildParentInterpolator[TreeLevel][child]);
}
}
/**
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations.
* This is a more optimal version : complexity p^4.
* This function handles extended cells.
*
*/
void initTensorM2MandL2L(const int TreeLevel, const FReal ParentWidth)
{
FReal ChildCoords[3][ORDER];
FPoint ChildCenter;
// 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
for (unsigned int child=0; child<8; ++child) {
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter, ExtendedCellRatio);
FChebTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory
ChildParentInterpolator[TreeLevel][child] = new FReal [3 * ORDER*ORDER];
assembleInterpolator(ORDER, ChildCoords[0], ChildParentInterpolator[TreeLevel][child]);
assembleInterpolator(ORDER, ChildCoords[1], ChildParentInterpolator[TreeLevel][child] + 1 * ORDER*ORDER);
assembleInterpolator(ORDER, ChildCoords[2], ChildParentInterpolator[TreeLevel][child] + 2 * ORDER*ORDER);
}
// init permutations
for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int index = k*ORDER*ORDER + j*ORDER + i;
perm[0][index] = k*ORDER*ORDER + j*ORDER + i;
perm[1][index] = i*ORDER*ORDER + k*ORDER + j;
perm[2][index] = j*ORDER*ORDER + i*ORDER + k;
}
}
// init permutations
for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int index = k*ORDER*ORDER + j*ORDER + i;
perm[0][index] = k*ORDER*ORDER + j*ORDER + i;
perm[1][index] = i*ORDER*ORDER + k*ORDER + j;
perm[2][index] = j*ORDER*ORDER + i*ORDER + k;
}
}
}
}
public:
/**
* Constructor: Initialize the Chebyshev polynomials at the Chebyshev
* 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 FChebInterpolator()
explicit FChebInterpolator(const int inTreeHeight=3,
const FReal inRootCellWidth=FReal(1.),
const FReal inCellWidthExtension=FReal(0.))
: TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth),
CellWidthExtension(inCellWidthExtension)
{
// initialize chebyshev polynomials of root nodes: T_o(x_j)
for (unsigned int o=1; o<ORDER; ++o)
......@@ -353,10 +383,35 @@ public:
// initialize root node ids
TensorType::setNodeIds(node_ids);
// initialize interpolation operator for non M2M and L2L (non leaf
// operations)
//this -> initM2MandL2L(); // non tensor-product interpolation
this -> initTensorM2MandL2L(); // tensor-product interpolation
// initialize interpolation operator for M2M/L2L (non leaf operations)
// 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
}
}
......@@ -365,8 +420,10 @@ public:
*/
~FChebInterpolator()
{
for (unsigned int l=0; l<TreeHeight; ++l)
for (unsigned int child=0; child<8; ++child)
delete [] ChildParentInterpolator[child];
if(ChildParentInterpolator[l][child] != nullptr)
delete [] ChildParentInterpolator[l][child];
}
......@@ -437,8 +494,6 @@ public:
const FReal *const * getChildParentInterpolator() const
{ return ChildParentInterpolator; }
const unsigned int * getPermutationsM2ML2L(unsigned int i) const
{ return perm[i]; }
......@@ -510,27 +565,27 @@ public:
*/
void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion,
FReal *const ParentExpansion) const
const FReal *const ChildExpansion,
FReal *const ParentExpansion,
const unsigned int TreeLevel = 2) const
{
FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER,
ChildParentInterpolator[TreeLevel][ChildIndex], ORDER,
const_cast<FReal*>(ChildExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-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);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-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);
for (unsigned int n=0; n<nnodes; ++n) ParentExpansion[perm[2][n]] += PermExp[n];
......@@ -538,25 +593,26 @@ public:
void applyL2L(const unsigned int ChildIndex,
const FReal *const ParentExpansion,
FReal *const ChildExpansion) const
const FReal *const ParentExpansion,
FReal *const ChildExpansion,
const unsigned int TreeLevel = 2) const
{
FReal Exp[nnodes], PermExp[nnodes];
// ORDER*ORDER*ORDER * (2*ORDER-1)
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER,
ChildParentInterpolator[TreeLevel][ChildIndex], ORDER,
const_cast<FReal*>(ParentExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
// ORDER*ORDER*ORDER * (2*ORDER-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);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
// ORDER*ORDER*ORDER * (2*ORDER-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);
for (unsigned int n=0; n<nnodes; ++n) ChildExpansion[perm[2][n]] += PermExp[n];
......
......@@ -71,12 +71,14 @@ public:
FChebTensorialKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal inBoxWidthExtension,
const FReal Epsilon,
const double inMatParam = 0.0)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inMatParam),
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight,
inBoxWidth,
inBoxWidthExtension,
Epsilon))
{ }
......@@ -85,11 +87,13 @@ public:
const ContainerClass* const SourceParticles)
{
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension);
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getMultipole(idxV*nRhs), SourceParticles);
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
......@@ -106,7 +110,7 @@ public:
void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/)
const int TreeLevel)
{
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
......@@ -117,8 +121,10 @@ public:
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxMul));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxMul),
ParentCell->getMultipole(idxMul));
AbstractBaseClass::Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(idxMul),
ParentCell->getMultipole(idxMul),
TreeLevel/*Cell width extension specific*/);
}
}
......@@ -134,8 +140,10 @@ public:
const int /*NumSourceCells*/,
const int TreeLevel)
{
// scale factor for homogeneous case
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 idxLhs=0; idxLhs < nLhs; ++idxLhs){
......@@ -167,7 +175,7 @@ public:
void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/)
const int TreeLevel)
{
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
......@@ -179,7 +187,10 @@ public:
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++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*/);
}
}
}//NLHS
......@@ -191,6 +202,8 @@ public:
ContainerClass* const TargetParticles)
{
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension);
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
......@@ -213,7 +226,7 @@ public:
// TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getLocal(idxV*nLhs), TargetParticles);
}
}
......
......@@ -66,15 +66,9 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
* size \f$\ell^3\times r\f$, and \f$316\f$ \f$C_t\f$, each of size \f$r\times
* r\f$.
*
* PB: FChebM2LHandler does not seem to support non_homogeneous kernels!
* In fact nothing appears to handle this here (i.e. adapt scaling and storage
* to MatrixKernelClass::Type). Given the relatively important cost of the
* Chebyshev variant, it is probably a choice not to have implemented this
* feature here but instead in the ChebyshevSym variant. But what if the
* kernel is non homogeneous and non symmetric (e.g. Dislocations)...
*
* TODO Specialize class (see UnifM2LHandler) OR prevent from using this
* class with non homogeneous kernels ?!
* PB: BEWARE! Homogeneous matrix kernels do not support cell width extension
* yet. Is it possible to find a reference width and a scale factor such that
* only 1 set of M2L ops can be used for all levels??
*
* @tparam ORDER interpolation order \f$\ell\f$
*/
......@@ -94,6 +88,8 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
FReal *U, *B;
FReal** C;
const FReal CellWidthExtension; //<! extension of cells width
const FReal epsilon; //<! accuracy which determines trucation of SVD
unsigned int rank; //<! truncation rank, satisfies @p epsilon
......@@ -109,12 +105,13 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
public:
FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal _epsilon)
: U(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal inCellWidthExtension, const FReal _epsilon)
: U(nullptr), B(nullptr), CellWidthExtension(inCellWidthExtension),
epsilon(_epsilon), rank(0)
{
// measure time
FTic time; time.tic();
// check if aready set
// check if already set
if (U||B) throw std::runtime_error("U or B operator already set");
// allocate C
......@@ -126,8 +123,11 @@ public:
if (C[d]) throw std::runtime_error("Compressed M2L operator already set");
// 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.);
rank = ComputeAndCompress<order>(MatrixKernel, ReferenceCellWidth, epsilon, U, C, B);
rank = ComputeAndCompress<order>(MatrixKernel, ReferenceCellWidth, 0., epsilon, U, C, B);
unsigned long sizeM2L = 343*ncmp*rank*rank*sizeof(FReal);
......@@ -210,8 +210,9 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
FReal **U, **B;
FReal*** C;
const unsigned int TreeHeight;
const FReal RootCellWidth;
const unsigned int TreeHeight; //<! number of levels
const FReal RootCellWidth; //<! width of root cell
const FReal CellWidthExtension; //<! extension of cells width
const FReal epsilon; //<! accuracy which determines trucation of SVD
unsigned int *rank; //<! truncation rank, satisfies @p epsilon
......@@ -228,9 +229,10 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
public:
FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal _epsilon)
FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal inCellWidthExtension, const FReal _epsilon)
: TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth),
CellWidthExtension(inCellWidthExtension),
epsilon(_epsilon)
{
// measure time
......@@ -265,7 +267,9 @@ public:
CellWidth /= FReal(2.); // at level 2
rank[0]=rank[1]=0;
for (unsigned int l=2; l<TreeHeight; ++l) {
rank[l] = ComputeAndCompress<order>(MatrixKernel, CellWidth, epsilon, U[l], C[l], B[l]);
// compute m2l operator on extended cell
rank[l] = ComputeAndCompress<order>(MatrixKernel, CellWidth, CellWidthExtension, epsilon, U[l], C[l], B[l]);
// update cell width
CellWidth /= FReal(2.); // at level l+1
}
unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*rank[2]*rank[2]*sizeof(FReal);
......@@ -354,6 +358,7 @@ public:
template <int ORDER, class MatrixKernelClass>
unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const FReal CellWidth,
const FReal CellWidthExtension,
const FReal epsilon,
FReal* &U,
FReal** &C,
......@@ -372,7 +377,8 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
// 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.), CellWidth, X);
const FReal ExtendedCellWidth(CellWidth+CellWidthExtension);
FChebTensor<order>::setRoots(FPoint(0.,0.,0.), ExtendedCellWidth, X);
// allocate memory and compute 316 m2l operators
FReal** _C;
......@@ -387,14 +393,13 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
// set roots of source cell (Y)
const FPoint cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
FChebTensor<order>::setRoots(cy, CellWidth, Y);
FChebTensor<order>::setRoots(cy, ExtendedCellWidth, Y);
// evaluate m2l operator
for (unsigned int n=0; n<nnodes; ++n)
for (unsigned int m=0; m<nnodes; ++m){
// _C[counter*nnodes*nnodes + n*nnodes + m]
// = MatrixKernel.evaluate(X[m], Y[n]);
// Compute current M2L interaction (block matrix)
// Compute current M2L interaction (block matrix)
FReal* block;
block = new FReal[ncmp];
MatrixKernel->evaluateBlock(X[m], Y[n], block);
......
......@@ -52,6 +52,10 @@ enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
* The table applyTab provides the indices in the reduced storage table
* corresponding to the application scheme depicted ealier.
*
* PB: BEWARE! Homogeneous matrix kernels do not support cell width extension
* yet. Is it possible to find a reference width and a scale factor such that
* only 1 set of M2L ops can be used for all levels??
*
*/
struct FInterpAbstractMatrixKernel : FNoCopyable
{
......
......@@ -128,10 +128,12 @@ public:
*
* @param[in] ChildIndex index of child according to Morton index
* @param[out] center
* @param[in] ExtendedCellRatio ratio between extended child and parent widths
*/
static
void setRelativeChildCenter(const unsigned int ChildIndex,
FPoint& ChildCenter)
FPoint& ChildCenter,
const FReal ExtendedCellRatio=FReal(.5))
{