Commit 56004d65 authored by Marc-Alexandre Espiaut's avatar Marc-Alexandre Espiaut
Browse files
parents f1be5e2d 134b371f
......@@ -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))