Commit 3ebd26e3 authored by Quentin Khan's avatar Quentin Khan

FChebCell/FUnifCell: simplify code

The multipole and local expansion share the same implementation. Their
types only differ by a template tag parameter. Consquently, the
`getMultipole`/`getLocal` methods are changed to `get`.

The method change is propageted to the files:

  - FChebKernel,
  - FChebSymKernel,
  - FChebTensorialKernel,
  - FAdaptChebKernel,
  - FUnifKernel,
  - FUnifTensorialKernel,
  - FAdaptUnifKernel,
  - utestChebyshevDirectPeriodic,
  - utestInterpolationMultiRhs,
parent fb301fa8
......@@ -50,7 +50,7 @@ public:
const FPoint<FReal> LeafCellCenter =
FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
FBase::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
LeafCell->getMultipole(0), SourceParticles);
LeafCell->get(0), SourceParticles);
}
template<class SymbolicData>
......@@ -63,7 +63,7 @@ public:
FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
// apply Sx and Px (grad Sx)
FBase::Interpolator->applyL2PTotal(LeafCellCenter, leafBoxWidth,
LeafCell->getLocal(0), TargetParticles);
LeafCell->get(0), TargetParticles);
}
template<class SymbolicData>
......@@ -117,7 +117,7 @@ public:
* physicalValues[idxPart];
}
local->getLocal(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);
local->get(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);
// Compute the last array elements one by one if they exist
if(idxPart < ((particles->getNbParticles() + FRealCount - 1) / FRealCount)) {
......@@ -129,7 +129,7 @@ public:
idxPart < static_cast<std::size_t>(particles->getNbParticles());
++idxPart)
{
local->getLocal(idxRhs)[m] +=
local->get(idxRhs)[m] +=
FBase::MatrixKernel->evaluate(
Xx, Xy, Xz,
pX[idxPart], pY[idxPart], pZ[idxPart])
......@@ -175,7 +175,7 @@ public:
for (unsigned int n=0; n<FBase::nnodes; ++n){
ComputeClass MultipoleExpansion =
FMath::ConvertTo<ComputeClass, FReal>(pole->getMultipole(idxRhs)[n]);
FMath::ConvertTo<ComputeClass, FReal>(pole->get(idxRhs)[n]);
ComputeClass YX = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getX());
ComputeClass YY = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getY());
......
......@@ -60,12 +60,12 @@ public:
// static_cast<int>(LeafSymbData->depth)));
// FBase::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
// LeafMultipole->getMultipole(0), SourceParticles);
// LeafMultipole->get(0), SourceParticles);
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // 2) apply Discrete Fourier Transform
// this->M2LHandler.applyZeroPaddingAndDFT(LeafMultipole->getMultipole(idxRhs),
// LeafMultipole->getTransformedMultipole(idxRhs));
// this->M2LHandler.applyZeroPaddingAndDFT(LeafMultipole->get(idxRhs),
// LeafMultipole->getTransformed(idxRhs));
// }
// }
......@@ -90,9 +90,9 @@ public:
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // 1) Apply Inverse Discete Fourier Transform
// FBase::M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs),
// FBase::M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformed(idxRhs),
// localExp + idxRhs*FBase::nnodes);
// FBlas::add(FBase::nnodes,const_cast<FReal*>(LeafCell->getLocal(idxRhs)),localExp + idxRhs*FBase::nnodes);
// FBlas::add(FBase::nnodes,const_cast<FReal*>(LeafCell->get(idxRhs)),localExp + idxRhs*FBase::nnodes);
// }
// // 2.a) apply Sx
// FBase::Interpolator->applyL2P(LeafCellCenter, leafBoxWidth,
......@@ -157,7 +157,7 @@ public:
* physicalValues[idxPart];
}
local->getLocal(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);
local->get(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);
// Compute the last array elements one by one if they exist
if(idxPart < ((particles->getNbParticles() + FRealCount - 1) / FRealCount)) {
......@@ -169,7 +169,7 @@ public:
idxPart < static_cast<std::size_t>(particles->getNbParticles());
++idxPart)
{
local->getLocal(idxRhs)[m] +=
local->get(idxRhs)[m] +=
FBase::MatrixKernel->evaluate(
Xx, Xy, Xz,
pX[idxPart], pY[idxPart], pZ[idxPart])
......@@ -221,7 +221,7 @@ public:
for (unsigned int n=0; n<FBase::nnodes; ++n){
ComputeClass MultipoleExpansion =
FMath::ConvertTo<ComputeClass, FReal>(pole->getMultipole(idxRhs)[n]);
FMath::ConvertTo<ComputeClass, FReal>(pole->get(idxRhs)[n]);
ComputeClass YX = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getX());
ComputeClass YY = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getY());
......
......@@ -39,13 +39,14 @@ class FChebCell : public FBasicCell, public FAbstractSendable
public:
struct multipole_t {
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
template<class Tag, std::size_t N>
struct exp_impl {
FReal exp[N * NVALS * VectorSize];
const FReal* getMultipole(const int inRhs) const
{ return this->multipole_exp + inRhs*VectorSize; }
FReal* getMultipole(const int inRhs)
{ return this->multipole_exp + inRhs*VectorSize; }
const FReal* get(const int inRhs) const
{ return this->exp + inRhs*VectorSize; }
FReal* get(const int inRhs)
{ return this->exp + inRhs*VectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
......@@ -54,45 +55,27 @@ public:
// to extend FAbstractSendable
template <class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const{
buffer.write(this->multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(this->exp, VectorSize*NVALS*NRHS);
}
template <class BufferReaderClass>
void deserialize(BufferReaderClass& buffer){
buffer.fillArray(this->multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(this->exp, VectorSize*NVALS*NRHS);
}
void reset() {
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(this->exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
}
};
struct local_expansion_t {
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
const FReal* getLocal(const int inRhs) const
{ return this->local_exp + inRhs*VectorSize; }
FReal* getLocal(const int inRhs)
{ return this->local_exp + inRhs*VectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
FSize getSavedSize() const {
return FSize(sizeof(FReal)) * VectorSize * N * NVALS;
}
// to extend FAbstractSendable
template <class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const{
buffer.write(this->local_exp, VectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void deserialize(BufferReaderClass& buffer){
buffer.fillArray(this->local_exp, VectorSize*NVALS*NLHS);
}
void reset() {
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
}
};
using multipole_t = exp_impl<class multipole_tag, NRHS>;
using local_expansion_t = exp_impl<class local_expansion_tag, NLHS>;
multipole_t m_data {};
local_expansion_t l_data {};
......@@ -187,7 +170,7 @@ public:
// const void print() const{
output <<" Multipole exp NRHS " <<NRHS <<" NVALS " <<NVALS << " VectorSize/2 " << cell.getVectorSize() *0.5<< std::endl;
for (int rhs= 0 ; rhs < NRHS ; ++rhs) {
const FReal* pole = cell.getMultipole(rhs);
const FReal* pole = cell.get(rhs);
for (int val= 0 ; val < NVALS ; ++val) {
output<< " val : " << val << " exp: " ;
for (int i= 0 ; i < cell.getVectorSize()/2 ; ++i) {
......
......@@ -97,12 +97,12 @@ public:
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
LeafMultipole->getMultipole(0), SourceParticles);
LeafMultipole->get(0), SourceParticles);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 2) apply B
M2LHandler->applyB(LeafMultipole->getMultipole(idxRhs),
LeafMultipole->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
M2LHandler->applyB(LeafMultipole->get(idxRhs),
LeafMultipole->get(idxRhs) + AbstractBaseClass::nnodes);
}
}
......@@ -118,13 +118,13 @@ public:
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildMultipoles[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(
ChildIndex, ChildMultipoles[ChildIndex]->getMultipole(idxRhs),
ParentMultipole->getMultipole(idxRhs));
ChildIndex, ChildMultipoles[ChildIndex]->get(idxRhs),
ParentMultipole->get(idxRhs));
}
}
// 2) apply B
M2LHandler->applyB(ParentMultipole->getMultipole(idxRhs),
ParentMultipole->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
M2LHandler->applyB(ParentMultipole->get(idxRhs),
ParentMultipole->get(idxRhs) + AbstractBaseClass::nnodes);
}
}
......@@ -142,8 +142,8 @@ public:
// cy.getY()-cx.getY(),
// cy.getZ()-cx.getZ()};
// M2LHandler->applyC(transfer, CellWidth,
// SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
// TargetCell->getLocal() + AbstractBaseClass::nnodes);
// SourceCells[idx]->get() + AbstractBaseClass::nnodes,
// TargetCell->get() + AbstractBaseClass::nnodes);
// }
// }
......@@ -156,12 +156,12 @@ public:
const int inSize)
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FReal *const CompressedLocalExpansion = TargetExpansion->getLocal(idxRhs) + AbstractBaseClass::nnodes;
FReal *const CompressedLocalExpansion = TargetExpansion->get(idxRhs) + AbstractBaseClass::nnodes;
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(1 << TargetSymb->getLevel()));
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idx = neighborPositions[idxExistingNeigh];
M2LHandler->applyC(idx, CellWidth,
SourceMultipoles[idxExistingNeigh]->getMultipole(idxRhs)
SourceMultipoles[idxExistingNeigh]->get(idxRhs)
+ AbstractBaseClass::nnodes,
CompressedLocalExpansion);
}
......@@ -175,10 +175,10 @@ public:
// const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
// for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
// const int idx = neighborPositions[idxExistingNeigh];
// FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idxExistingNeigh]->getMultipole())+AbstractBaseClass::nnodes,
// FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idxExistingNeigh]->get())+AbstractBaseClass::nnodes,
// MultipoleExpansion+idx*rank);
//
// M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
// M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->get() + AbstractBaseClass::nnodes);
// }
......@@ -190,15 +190,15 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
M2LHandler->applyU(ParentExpansion->getLocal(idxRhs) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(ParentExpansion)->getLocal(idxRhs));
M2LHandler->applyU(ParentExpansion->get(idxRhs) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(ParentExpansion)->get(idxRhs));
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildExpansions[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(
ChildIndex,
ParentExpansion->getLocal(idxRhs),
ChildExpansions[ChildIndex]->getLocal(idxRhs));
ParentExpansion->get(idxRhs),
ChildExpansions[ChildIndex]->get(idxRhs));
}
}
}
......@@ -218,24 +218,24 @@ public:
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
M2LHandler->applyU(LeafExpansion->getLocal(idxRhs) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(LeafExpansion)->getLocal(idxRhs));
M2LHandler->applyU(LeafExpansion->get(idxRhs) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(LeafExpansion)->get(idxRhs));
}
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(0),
// LeafCell->get(0),
// TargetParticles);
//// 2.b) apply Px (grad Sx)
//AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(0),
// LeafCell->get(0),
// TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, leafBoxWidth,
LeafExpansion->getLocal(0), TargetParticles);
LeafExpansion->get(0), TargetParticles);
}
......
......@@ -191,7 +191,7 @@ public:
FReal leafBoxWidth = AbstractBaseClass::BoxWidth / FReal(1 << LeafSymbData->getLevel());
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
LeafCell->getMultipole(0), SourceParticles);
LeafCell->get(0), SourceParticles);
}
......@@ -204,13 +204,13 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// Reset the Parent expansion to zero
// FBlas::scal(nnodes*2, FReal(0.), ParentMultipole->getMultipole(idxRhs));
// FBlas::scal(nnodes*2, FReal(0.), ParentMultipole->get(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
// apply Sy
if (ChildMultipoles[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(
ChildIndex, ChildMultipoles[ChildIndex]->getMultipole(idxRhs),
ParentMultipole->getMultipole(idxRhs));
ChildIndex, ChildMultipoles[ChildIndex]->get(idxRhs),
ParentMultipole->get(idxRhs));
}
}
}
......@@ -238,7 +238,7 @@ public:
const unsigned int count = (countExp[pidx])++;
FReal *const mul = Mul[pidx] + count*nnodes;
const unsigned int *const pvec = SymHandler->pvectors[idx];
const FReal *const MultiExp = SourceMultipoles[idxExistingNeigh]->getMultipole(idxRhs);
const FReal *const MultiExp = SourceMultipoles[idxExistingNeigh]->get(idxRhs);
/*
// no loop unrolling
......@@ -311,7 +311,7 @@ public:
#endif
// permute and add contribution to local expansions
FReal *const LocalExpansion = TargetExpansion->getLocal(idxRhs);
FReal *const LocalExpansion = TargetExpansion->get(idxRhs);
memset(countExp, 0, sizeof(int) * 343);
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idx = neighborPositions[idxExistingNeigh];
......@@ -356,7 +356,7 @@ public:
const int idx = neighborPositions[idxExistingNeigh];
const unsigned int pidx = SymHandler->pindices[idx];
const unsigned int count = (countExp[pidx])++;
const FReal *const MultiExp = SourceCells[idxExistingNeigh]->getMultipole();
const FReal *const MultiExp = SourceCells[idxExistingNeigh]->get();
for (unsigned int n=0; n<nnodes; ++n)
Mul[pidx][count*nnodes + SymHandler->pvectors[idx][n]] = MultiExp[n];
}
......@@ -377,7 +377,7 @@ public:
}
// permute and add contribution to local expansions
FReal *const LocalExpansion = TargetCell->getLocal();
FReal *const LocalExpansion = TargetCell->get();
memset(countExp, 0, sizeof(int) * 343);
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idx = neighborPositions[idxExistingNeigh];
......@@ -393,7 +393,7 @@ public:
/*
void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
const int neighborPositions[], const int inSize, const int TreeLevel) override {
FReal *const LocalExpansion = TargetCell->getLocal();
FReal *const LocalExpansion = TargetCell->get();
const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
......@@ -407,7 +407,7 @@ public:
const unsigned int idx = ((i+3) * 7 + (j+3)) * 7 + (k+3);
if (SourceCells[idx]) {
const FReal *const MultiExp = SourceCells[idx]->getMultipole();
const FReal *const MultiExp = SourceCells[idx]->get();
// permute
for (unsigned int n=0; n<nnodes; ++n) PermMultiExp[pvectors[idx][n]] = MultiExp[n];
......@@ -441,8 +441,8 @@ public:
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildExpansions[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(
ChildIndex, ParentExpansion->getLocal(idxRhs),
ChildExpansions[ChildIndex]->getLocal(idxRhs));
ChildIndex, ParentExpansion->get(idxRhs),
ChildExpansions[ChildIndex]->get(idxRhs));
}
}
}
......@@ -458,17 +458,17 @@ public:
// // a) apply Sx
// AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// LeafCell->get(),
// TargetParticles);
// // b) apply Px (grad Sx)
// AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// LeafCell->get(),
// TargetParticles);
// c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(0), TargetParticles);
LeafCell->get(0), TargetParticles);
}
......
......@@ -102,7 +102,7 @@ public:
AbstractBaseClass::Interpolator->applyP2M(
LeafCellCenter,
ExtendedLeafCellWidth,
LeafMultipole->getMultipole(idxV*nRhs),
LeafMultipole->get(idxV*nRhs),
SourceParticles);
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
......@@ -110,8 +110,8 @@ public:
int idxMul = idxV*nRhs + idxRhs;
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler->applyB(LeafMultipole->getMultipole(idxMul),
LeafMultipole->getMultipole(idxMul) + AbstractBaseClass::nnodes);
M2LHandler->applyB(LeafMultipole->get(idxMul),
LeafMultipole->get(idxMul) + AbstractBaseClass::nnodes);
}
}
......@@ -131,20 +131,20 @@ public:
int idxMul = idxV*nRhs + idxRhs;
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentMultipole->getMultipole(idxMul));
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentMultipole->get(idxMul));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildMultipoles[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(
ChildIndex,
ChildMultipoles[ChildIndex]->getMultipole(idxMul),
ParentMultipole->getMultipole(idxMul),
ChildMultipoles[ChildIndex]->get(idxMul),
ParentMultipole->get(idxMul),
TreeLevel/*Cell width extension specific*/);
}
}
// 2) apply B (PB: Tensorial version is just a basic copy)
M2LHandler->applyB(ParentMultipole->getMultipole(idxMul),
ParentMultipole->getMultipole(idxMul) + AbstractBaseClass::nnodes);
M2LHandler->applyB(ParentMultipole->get(idxMul),
ParentMultipole->get(idxMul) + AbstractBaseClass::nnodes);
}
}
......@@ -169,7 +169,7 @@ public:
// update local index
const int idxLoc = idxV*nLhs + idxLhs;
FReal *const CompressedLocalExpansion = TargetExpansion->getLocal(idxLoc) + AbstractBaseClass::nnodes;
FReal *const CompressedLocalExpansion = TargetExpansion->get(idxLoc) + AbstractBaseClass::nnodes;
// update idxRhs
const int idxRhs = idxLhs % nPv;
......@@ -183,7 +183,7 @@ public:
const int idx = neighborPositions[idxExistingNeigh];
M2LHandler->applyC(
idx, TreeLevel, scale, d,
SourceMultipoles[idxExistingNeigh]->getMultipole(idxMul) + AbstractBaseClass::nnodes,
SourceMultipoles[idxExistingNeigh]->get(idxMul) + AbstractBaseClass::nnodes,
CompressedLocalExpansion);
}
}// NLHS=NPOT*NPV
......@@ -204,14 +204,14 @@ public:
// 1) apply U (PB: Tensorial version is just a basic copy)
M2LHandler->applyU(
ParentExpansion->getLocal(idxLoc) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(ParentExpansion)->getLocal(idxLoc));
ParentExpansion->get(idxLoc) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(ParentExpansion)->get(idxLoc));
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildExpansions[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(ChildIndex,
ParentExpansion->getLocal(idxLoc),
ChildExpansions[ChildIndex]->getLocal(idxLoc),
ParentExpansion->get(idxLoc),
ChildExpansions[ChildIndex]->get(idxLoc),
TreeLevel/*Cell width extension specific*/);
}
}
......@@ -234,13 +234,13 @@ public:
int idxLoc = idxV*nLhs + idxLhs;
// 1) apply U (PB: Tensorial version is just a basic copy)
M2LHandler->applyU(
LeafExpansion->getLocal(idxLoc) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(LeafExpansion)->getLocal(idxLoc));
LeafExpansion->get(idxLoc) + AbstractBaseClass::nnodes,
const_cast<typename CellClass::local_expansion_t*>(LeafExpansion)->get(idxLoc));
}
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, ExtendedLeafCellWidth,
LeafExpansion->getLocal(idxV*nLhs), TargetParticles);
LeafExpansion->get(idxV*nLhs), TargetParticles);
}
}
......
......@@ -44,19 +44,21 @@ class FUnifCell : public FBasicCell, public FAbstractSendable
static const int TransformedVectorSize = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1);
public:
struct multipole_t {
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
template<class Tag, std::size_t N>
struct exp_impl {
FReal exp[N * NVALS * VectorSize]; //< Multipole expansion
// Multipole expansion in Fourier space
FComplex<FReal> transformed_multipole_exp[NRHS * NVALS * TransformedVectorSize];
FComplex<FReal> transformed_exp[N * NVALS * TransformedVectorSize];
const FReal* getMultipole(const int inRhs) const
{ return this->multipole_exp + inRhs*VectorSize; }
FReal* getMultipole(const int inRhs)
{ return this->multipole_exp + inRhs*VectorSize; }
const FComplex<FReal>* getTransformedMultipole(const int inRhs) const
{ return this->transformed_multipole_exp + inRhs*TransformedVectorSize; }
FComplex<FReal>* getTransformedMultipole(const int inRhs)
{ return this->transformed_multipole_exp + inRhs*TransformedVectorSize; }
const FReal* get(const int inRhs) const
{ return this->exp + inRhs*VectorSize; }
FReal* get(const int inRhs)
{ return this->exp + inRhs*VectorSize; }
const FComplex<FReal>* getTransformed(const int inRhs) const
{ return this->transformed_exp + inRhs*TransformedVectorSize; }
FComplex<FReal>* getTransformed(const int inRhs)
{ return this->transformed_exp + inRhs*TransformedVectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
......@@ -65,89 +67,63 @@ public:
// to extend FAbstractSendable
template <class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const{
buffer.write(this->multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(this->transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
buffer.write(this->exp, VectorSize*NVALS*N);
buffer.write(this->transformed_exp, TransformedVectorSize*NVALS*N);
}
template <class BufferReaderClass>
void deserialize(BufferReaderClass& buffer){
buffer.fillArray(this->multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(this->transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
}
};
struct local_expansion_t {
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
// Local expansion in Fourier space
FComplex<FReal> transformed_local_exp[NLHS * NVALS * TransformedVectorSize];
const FReal* getLocal(const int inRhs) const
{ return this->local_exp + inRhs*VectorSize; }
FReal* getLocal(const int inRhs)
{ return this->local_exp + inRhs*VectorSize; }
const FComplex<FReal>* getTransformedLocal(const int inRhs) const
{ return this->transformed_local_exp + inRhs*TransformedVectorSize; }
FComplex<FReal>* getTransformedLocal(const int inRhs)
{ return this->transformed_local_exp + inRhs*TransformedVectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
buffer.fillArray(this->exp, VectorSize*NVALS*N);
buffer.fillArray(this->transformed_exp, TransformedVectorSize*NVALS*N);
}
// to extend FAbstractSendable
///////////////////////////////////////////////////////
// to extend Serializable
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const{
buffer.write(this->local_exp, VectorSize*NVALS*NLHS);
buffer.write(this->transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
void save(BufferWriterClass& buffer) const{
buffer.write(this->exp, VectorSize*NVALS*NRHS);
buffer.write(this->transformed_exp, TransformedVectorSize*NVALS*NRHS);
}
template <class BufferReaderClass>
void deserialize(BufferReaderClass& buffer){
buffer.fillArray(this->local_exp, VectorSize*NVALS*NLHS);
buffer.fillArray(this->transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
void restore(BufferReaderClass& buffer){
buffer.fillArray(this->exp, VectorSize*NVALS*NRHS);
buffer.fillArray(this->transformed_exp, TransformedVectorSize*NVALS*NRHS);
}
FSize getSavedSize() const {
return N * NVALS * VectorSize * (FSize) sizeof(FReal)
+ N * NVALS * TransformedVectorSize * (FSize) sizeof(FComplex<FReal>);
}
};
using multipole_t = exp_impl<class multipole_tag, NRHS>;
using local_expansion_t = exp_impl<class local_expansion_tag, NLHS>;
multipole_t m_data {};
local_expansion_t l_data {};
multipole_t* up = &m_data;
local_expansion_t* down = &l_data;
FUnifCell(){
memset(up->multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(down->local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
memset(up->transformed_multipole_exp, 0,
sizeof(FComplex<FReal>) * NRHS * NVALS * TransformedVectorSize);
memset(down->transformed_local_exp, 0,
sizeof(FComplex<FReal>) * NLHS * NVALS * TransformedVectorSize);
}
~FUnifCell() {}
bool hasMultipoleData() const noexcept {
return up != nullptr;
return true;
}
bool hasLocalExpansionData() const noexcept {
return up != nullptr;
return true;
}
multipole_t& getMultipoleData() noexcept {
return *up;
return m_data;
}
const multipole_t& getMultipoleData() const noexcept {
return *up;
return m_data;
}
local_expansion_t& getLocalExpansionData() noexcept {
return *down;
return l_data;
}
const local_expansion_t& getLocalExpansionData() const noexcept {
return *down;
return l_data;
}
......@@ -156,22 +132,22 @@ public:
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const{
up->serialize(buffer);
m_data.serialize(buffer);
}
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
up->deserialize(buffer);
m_data.deserialize(buffer);
}
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const{
down->serialize(buffer);
l_data.serialize(buffer);
}
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
down->deserialize(buffer);
l_data.deserialize(buffer);
}
///////////////////////////////////////////////////////
......@@ -180,37 +156,36 @@ public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
buffer.write(up->multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(up->transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
buffer.write(down->local_exp, VectorSize*NVALS*NLHS);
buffer.write(down->transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
m_data.save(buffer);
l_data.save(buffer);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
buffer.fillArray(up->multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(up->transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
buffer.fillArray(down->local_exp, VectorSize*NVALS*NLHS);
buffer.fillArray(down->transformed_local_exp, TransformedVectorSize*NVALS*NLHS);