Skip to content
GitLab
Commit 95626445 authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files

Remove old fashion periodic system from spherical kernel

parent 6f55fa9a
Hide whitespace changes
Inline Side-by-side
Src/Kernels/Spherical/FAbstractSphericalKernel.hpp
View file @ 95626445
......@@ -33,8 +33,6 @@ protected:
const FReal boxWidth; //< the box width at leaf level
const int treeHeight; //< The height of the tree
const int periodicLevels; //< The number of levels above 1 used for periodicity
const FReal widthAtLeafLevel; //< the width of a box at leaf level
const FReal widthAtLeafLevelDiv2; //< the width of a box at leaf level divided by 2
const FPoint boxCorner; //< the corner of the box system
......@@ -48,15 +46,15 @@ protected:
/** Alloc and init pre-vectors*/
void allocAndInit(){
preL2LTransitions = new FComplexe*[treeHeight + periodicLevels];
memset(preL2LTransitions.getPtr(), 0, (treeHeight + periodicLevels) * sizeof(FComplexe*));
preM2MTransitions = new FComplexe*[treeHeight + periodicLevels];
memset(preM2MTransitions.getPtr(), 0, (treeHeight + periodicLevels) * sizeof(FComplexe*));
preL2LTransitions = new FComplexe*[treeHeight ];
memset(preL2LTransitions.getPtr(), 0, (treeHeight) * sizeof(FComplexe*));
preM2MTransitions = new FComplexe*[treeHeight];
memset(preM2MTransitions.getPtr(), 0, (treeHeight) * sizeof(FComplexe*));
FReal treeWidthAtLevel = (boxWidth * FReal(1 << periodicLevels))/2;
for(int idxLevel = -periodicLevels ; idxLevel < treeHeight - 1 ; ++idxLevel ){
preL2LTransitions[idxLevel + periodicLevels] = new FComplexe[ 8 * harmonic.getExpSize()];
preM2MTransitions[idxLevel + periodicLevels] = new FComplexe[ 8 * harmonic.getExpSize()];
FReal treeWidthAtLevel = (boxWidth)/2;
for(int idxLevel = 0 ; idxLevel < treeHeight - 1 ; ++idxLevel ){
preL2LTransitions[idxLevel] = new FComplexe[ 8 * harmonic.getExpSize()];
preM2MTransitions[idxLevel] = new FComplexe[ 8 * harmonic.getExpSize()];
const FPoint father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel);
treeWidthAtLevel /= 2;
......@@ -72,7 +70,7 @@ protected:
);
harmonic.computeInner(FSpherical(M2MVector));
FMemUtils::copyall<FComplexe>(&preM2MTransitions[idxLevel + periodicLevels][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
FMemUtils::copyall<FComplexe>(&preM2MTransitions[idxLevel][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
const FPoint L2LVector (
(treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))) - father.getX(),
......@@ -81,7 +79,7 @@ protected:
);
harmonic.computeInner(FSpherical(L2LVector));
FMemUtils::copyall<FComplexe>(&preL2LTransitions[idxLevel + periodicLevels][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
FMemUtils::copyall<FComplexe>(&preL2LTransitions[idxLevel][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
}
}
}
......@@ -97,11 +95,10 @@ protected:
public:
/** Kernel constructor */
FAbstractSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inPeriodicLevel = 0)
FAbstractSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
: devP(inDevP),
boxWidth(inBoxWidth),
treeHeight(inTreeHeight),
periodicLevels(inPeriodicLevel),
widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))),
widthAtLeafLevelDiv2(widthAtLeafLevel/2),
boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),inBoxCenter.getZ()-(inBoxWidth/2)),
......@@ -117,7 +114,6 @@ public:
: devP(other.devP),
boxWidth(other.boxWidth),
treeHeight(other.treeHeight),
periodicLevels(other.periodicLevels),
widthAtLeafLevel(other.widthAtLeafLevel),
widthAtLeafLevelDiv2(other.widthAtLeafLevelDiv2),
boxCorner(other.boxCorner),
......@@ -130,10 +126,10 @@ public:
/** Default destructor */
virtual ~FAbstractSphericalKernel(){
if(preL2LTransitions.isLast()){
FMemUtils::DeleteAll(preL2LTransitions.getPtr(), treeHeight + periodicLevels);
FMemUtils::DeleteAll(preL2LTransitions.getPtr(), treeHeight);
}
if(preM2MTransitions.isLast()){
FMemUtils::DeleteAll(preM2MTransitions.getPtr(), treeHeight + periodicLevels);
FMemUtils::DeleteAll(preM2MTransitions.getPtr(), treeHeight);
}
}
......@@ -155,7 +151,7 @@ public:
void M2M(CellClass* const FRestrict inPole, const CellClass *const FRestrict *const FRestrict inChild, const int inLevel) {
FComplexe* FRestrict const multipole_exp_target = inPole->getMultipole();
// iter on each child and process M2M
const FComplexe* FRestrict const preM2MTransitionsAtLevel = preM2MTransitions[inLevel + periodicLevels];
const FComplexe* FRestrict const preM2MTransitionsAtLevel = preM2MTransitions[inLevel];
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if(inChild[idxChild]){
multipoleToMultipole(multipole_exp_target, inChild[idxChild]->getMultipole(), &preM2MTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
......@@ -175,7 +171,7 @@ public:
}
std::cout <<std::endl;
//
const FComplexe* FRestrict const preM2MTransitionsAtLevel = preM2MTransitions[inLevel + periodicLevels];
const FComplexe* FRestrict const preM2MTransitionsAtLevel = preM2MTransitions[inLevel];
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if(inChild[idxChild]){
multipoleToMultipole(multipole_exp_target, inChild[idxChild]->getMultipole(), &preM2MTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
......@@ -195,7 +191,7 @@ public:
/** L2L with a cell and all its child */
void L2L(const CellClass* const FRestrict pole, CellClass* FRestrict *const FRestrict child, const int inLevel) {
// iter on each child and process L2L
const FComplexe* FRestrict const preL2LTransitionsAtLevel = preL2LTransitions[inLevel + periodicLevels];
const FComplexe* FRestrict const preL2LTransitionsAtLevel = preL2LTransitions[inLevel];
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if(child[idxChild]){
localToLocal(child[idxChild]->getLocal(), pole->getLocal(), &preL2LTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
......@@ -229,63 +225,6 @@ public:
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
ContainerClass* const directNeighborsParticles[27], const int /*size*/){
if( periodicLevels == 0 ){
P2PNoPeriodic(targets, sources, directNeighborsParticles);
}
else {
P2PPeriodic(inLeafPosition, targets, sources, directNeighborsParticles);
}
}
/** This P2P has to be used when target != sources
* It will proceed an direct interation no mutual
*
* It takes all the target particles from the current leaf,
* then it computes the sources/targets interaction in this leaf,
* then it computes the sources/targets inteactions between this leaf and the
* neighbors.
*/
void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,
ContainerClass* const directNeighborsParticles[27], const int /*size*/){
// TODO manage for periodic
for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < 27 ; ++idxDirectNeighbors){
if( directNeighborsParticles[idxDirectNeighbors] ){
// For all particles in current leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
ParticleClass target( iterTarget.data() );
// For all the particles in the other leaf
typename ContainerClass::BasicIterator iterSource(*directNeighborsParticles[idxDirectNeighbors]);
while( iterSource.hasNotFinished() ){
directInteraction(&target, iterSource.data());
iterSource.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
}
}
private:
///////////////////////////////////////////////////////////////////////////////
// P2P possibilities
///////////////////////////////////////////////////////////////////////////////
/** This P2P has to be used when target != sources
* It will proceed an direct interation no mutual
*
* It takes all the target particles from the current leaf,
* then it computes the sources/targets interaction in this leaf,
* then it computes the sources/targets inteactions between this leaf and the
* neighbors.
*/
void P2PNoPeriodic(ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
ContainerClass* const directNeighborsParticles[27]){
const bool isNotTsm = (targets != sources);
if( isNotTsm ) { // Compute interaction in this leaf
{
......@@ -366,7 +305,6 @@ private:
}
}
/** This P2P has to be used when target != sources
* It will proceed an direct interation no mutual
*
......@@ -375,113 +313,33 @@ private:
* then it computes the sources/targets inteactions between this leaf and the
* neighbors.
*/
void P2PPeriodic(const FTreeCoordinate& inLeafPosition,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
ContainerClass* const directNeighborsParticles[27]){
const int limite = FMath::pow2(treeHeight-1) - 1;
if( inLeafPosition.getX() == 0 || inLeafPosition.getX() == limite ||
inLeafPosition.getY() == 0 || inLeafPosition.getY() == limite ||
inLeafPosition.getZ() == 0 || inLeafPosition.getZ() == limite){
const bool isNotTsm = (targets != sources);
if( isNotTsm ) { // Compute interaction in this leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
// We copy the target particle to work with a particle in the heap
ParticleClass target( iterTarget.data() );
// For all the source particles in the same leaf
typename ContainerClass::ConstBasicIterator iterSameBox(*sources);
while( iterSameBox.hasNotFinished() ){
directInteraction(&target, iterSameBox.data());
iterSameBox.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
else { // Compute interaction in this leaf
void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,
ContainerClass* const directNeighborsParticles[27], const int /*size*/){
// TODO manage for periodic
for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < 27 ; ++idxDirectNeighbors){
if( directNeighborsParticles[idxDirectNeighbors] ){
// For all particles in current leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
// We copy the target particle to work with a particle in the heap
ParticleClass target( iterTarget.data() );
// For all the particles in the other leaf
typename ContainerClass::BasicIterator iterSource(*directNeighborsParticles[idxDirectNeighbors]);
while( iterSource.hasNotFinished() ){
directInteraction(&target, iterSource.data());
// For all particles after the current one
typename ContainerClass::BasicIterator iterSameBox = iterTarget;
iterSameBox.gotoNext();
while( iterSameBox.hasNotFinished() ){
directInteractionMutual(&target, &iterSameBox.data());
iterSameBox.gotoNext();
iterSource.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
{ // Compute interactions with other leaves
const int EndLoop = ( targets != sources ? 27 : 14);
// For all the neigbors leaves
for(int idxDirectNeighbors = 0 ; idxDirectNeighbors < EndLoop ; ++idxDirectNeighbors){
if( directNeighborsParticles[idxDirectNeighbors] ){
const int zrelatif = (idxDirectNeighbors%3) - 1;
const int yrelatif = ((idxDirectNeighbors%(3*3))/3)-1;
const int xrelatif = (idxDirectNeighbors / (3*3))-1;
FReal xoffset = 0;
FReal yoffset = 0;
FReal zoffset = 0;
if( inLeafPosition.getX() == 0 && xrelatif == -1){
xoffset = -boxWidth;
}
else if (inLeafPosition.getX() == limite && xrelatif == 1){
xoffset = boxWidth;
}
if( inLeafPosition.getY() == 0 && yrelatif == -1){
yoffset = -boxWidth;
}
else if (inLeafPosition.getY() == limite && yrelatif == 1){
yoffset = boxWidth;
}
if( inLeafPosition.getZ() == 0 && zrelatif == -1){
zoffset = -boxWidth;
}
else if (inLeafPosition.getZ() == limite && zrelatif == 1){
zoffset = boxWidth;
}
// For all particles in current leaf
typename ContainerClass::BasicIterator iterTarget(*targets);
while( iterTarget.hasNotFinished() ){
ParticleClass target( iterTarget.data() );
// For all the particles in the other leaf
typename ContainerClass::ConstBasicIterator iterSource(*directNeighborsParticles[idxDirectNeighbors]);
while( iterSource.hasNotFinished() ){
ParticleClass sourcePart( iterSource.data() );
sourcePart.setPosition(
sourcePart.getPosition().getX() + xoffset,
sourcePart.getPosition().getY() + yoffset,
sourcePart.getPosition().getZ() + zoffset
);
if(isNotTsm) directInteraction(&target, sourcePart);
else directInteractionMutual(&target, &sourcePart);;
iterSource.gotoNext();
}
// Set data and progress
iterTarget.setData(target);
iterTarget.gotoNext();
}
}
}
}
}
else{
P2PNoPeriodic(targets, sources, directNeighborsParticles);
}
}
private:
///////////////////////////////////////////////////////////////////////////////
// Computation
......
Src/Kernels/Spherical/FSphericalBlasKernel.hpp
View file @ 95626445
......@@ -47,13 +47,13 @@ protected:
// M2L transfer, there is a maximum of 3 neighbors in each direction,
// so 6 in each dimension
FReal treeWidthAtLevel = Parent::boxWidth * FReal( 1 << Parent::periodicLevels);
preM2LTransitions = new FComplexe**[Parent::treeHeight + Parent::periodicLevels];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe**) * (Parent::treeHeight + Parent::periodicLevels));
FReal treeWidthAtLevel = Parent::boxWidth;
preM2LTransitions = new FComplexe**[Parent::treeHeight];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe**) * (Parent::treeHeight));
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
preM2LTransitions[idxLevel + Parent::periodicLevels] = new FComplexe*[(7 * 7 * 7)];
memset(preM2LTransitions[idxLevel + Parent::periodicLevels], 0, sizeof(FComplexe*) * (7*7*7));
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
preM2LTransitions[idxLevel] = new FComplexe*[(7 * 7 * 7)];
memset(preM2LTransitions[idxLevel], 0, sizeof(FComplexe*) * (7*7*7));
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
......@@ -93,7 +93,7 @@ protected:
matrix[idxCol * FF_MATRIX_ROW_DIM + idxRow] = workMatrix[idxCol + idxRow * FF_MATRIX_COLUMN_DIM];
}
}
preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)] = matrix;
preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)] = matrix;
}
}
}
......@@ -113,8 +113,8 @@ public:
* @param inBoxWidth the size of the simulation box
* @param inPeriodicLevel the number of level upper to 0 that will be requiried
*/
FSphericalBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inPeriodicLevel = 0)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter, inPeriodicLevel),
FSphericalBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
FF_MATRIX_ROW_DIM(Parent::harmonic.getExpSize()), FF_MATRIX_COLUMN_DIM(Parent::harmonic.getNExpSize()),
FF_MATRIX_SIZE(FF_MATRIX_ROW_DIM * FF_MATRIX_COLUMN_DIM),
temporaryMultiSource(new FComplexe[FF_MATRIX_COLUMN_DIM]),
......@@ -136,16 +136,16 @@ public:
~FSphericalBlasKernel(){
delete[] temporaryMultiSource;
if(preM2LTransitions.isLast()){
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
delete[] preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)];
delete[] preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)];
}
}
}
}
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight + Parent::periodicLevels);
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight);
}
}
......@@ -155,7 +155,7 @@ public:
// For all neighbors compute M2L
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
if( distantNeighbors[idxNeigh] ){
const FComplexe* const transitionVector = preM2LTransitions[inLevel + Parent::periodicLevels][idxNeigh];
const FComplexe* const transitionVector = preM2LTransitions[inLevel][idxNeigh];
multipoleToLocal(pole->getLocal(), distantNeighbors[idxNeigh]->getMultipole(), transitionVector);
}
}
......
Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp
View file @ 95626445
......@@ -62,13 +62,13 @@ protected:
// M2L transfer, there is a maximum of 3 neighbors in each direction,
// so 6 in each dimension
FReal treeWidthAtLevel = Parent::boxWidth * FReal( 1 << Parent::periodicLevels);
preM2LTransitions = new FComplexe**[Parent::treeHeight + Parent::periodicLevels];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe**) * (Parent::treeHeight + Parent::periodicLevels));
FReal treeWidthAtLevel = Parent::boxWidth;
preM2LTransitions = new FComplexe**[Parent::treeHeight];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe**) * (Parent::treeHeight));
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
preM2LTransitions[idxLevel + Parent::periodicLevels] = new FComplexe*[(7 * 7 * 7)];
memset(preM2LTransitions[idxLevel + Parent::periodicLevels], 0, sizeof(FComplexe*) * (7*7*7));
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
preM2LTransitions[idxLevel] = new FComplexe*[(7 * 7 * 7)];
memset(preM2LTransitions[idxLevel], 0, sizeof(FComplexe*) * (7*7*7));
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
......@@ -108,7 +108,7 @@ protected:
matrix[idxCol * FF_MATRIX_ROW_DIM + idxRow] = workMatrix[idxCol + idxRow * FF_MATRIX_COLUMN_DIM];
}
}
preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)] = matrix;
preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)] = matrix;
}
}
}
......@@ -128,8 +128,8 @@ public:
* @param inBoxWidth the size of the simulation box
* @param inPeriodicLevel the number of level upper to 0 that will be requiried
*/
FSphericalBlockBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inBlockSize = 512, const int inPeriodicLevel = 0)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter, inPeriodicLevel),
FSphericalBlockBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inBlockSize = 512)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
FF_MATRIX_ROW_DIM(Parent::harmonic.getExpSize()), FF_MATRIX_COLUMN_DIM(Parent::harmonic.getNExpSize()),
FF_MATRIX_SIZE(FF_MATRIX_ROW_DIM * FF_MATRIX_COLUMN_DIM),
BlockSize(inBlockSize),
......@@ -156,16 +156,16 @@ public:
delete[] multipoleMatrix;
delete[] localMatrix;
if(preM2LTransitions.isLast()){
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
delete[] preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)];
delete[] preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)];
}
}
}
}
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight + Parent::periodicLevels);
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight);
}
}
......@@ -235,7 +235,7 @@ public:
*
*Remark: here we have always j+n >= |-k-l|
*
* const FComplexe*const M2LTransfer = preM2LTransitions[inLevel + Parent::periodicLevels][interactionIndex];
* const FComplexe*const M2LTransfer = preM2LTransitions[inLevel][interactionIndex];
* for(int idxK = 0 ; idxK < interactions[interactionIndex].getSize() ; ++idxK){
* for(int idxRow = 0 ; idxRow < FF_MATRIX_ROW_DIM ; ++idxRow){
* FComplexe compute;
......@@ -263,7 +263,7 @@ public:
FF_MATRIX_COLUMN_DIM,
interactions[interactionIndex].getSize(),
one,
FComplexe::ToFReal(preM2LTransitions[inLevel + Parent::periodicLevels][interactionIndex]),
FComplexe::ToFReal(preM2LTransitions[inLevel][interactionIndex]),
FF_MATRIX_ROW_DIM,
FComplexe::ToFReal(multipoleMatrix),
FF_MATRIX_COLUMN_DIM,
......
Src/Kernels/Spherical/FSphericalKernel.hpp
View file @ 95626445
......@@ -36,13 +36,13 @@ protected:
void allocAndInit(){
// M2L transfer, there is a maximum of 3 neighbors in each direction,
// so 6 in each dimension
preM2LTransitions = new FComplexe*[Parent::treeHeight + Parent::periodicLevels];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe*) * (Parent::treeHeight + Parent::periodicLevels));
preM2LTransitions = new FComplexe*[Parent::treeHeight];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe*) * (Parent::treeHeight));
// We start from the higher level
FReal treeWidthAtLevel = Parent::boxWidth * FReal(1 << Parent::periodicLevels);
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
FReal treeWidthAtLevel = Parent::boxWidth;
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
// Allocate data for this level
preM2LTransitions[idxLevel + Parent::periodicLevels] = new FComplexe[(7 * 7 * 7) * devM2lP];
preM2LTransitions[idxLevel] = new FComplexe[(7 * 7 * 7) * devM2lP];
// Precompute transfer vector
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
......@@ -50,7 +50,7 @@ protected:
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const FPoint relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
Parent::harmonic.computeOuter(FSpherical(relativePos));
FMemUtils::copyall<FComplexe>(&preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)], Parent::harmonic.result(), Parent::harmonic.getExpSize());
FMemUtils::copyall<FComplexe>(&preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)], Parent::harmonic.result(), Parent::harmonic.getExpSize());
}
}
}
......@@ -68,8 +68,8 @@ public:
* @param inBoxWidth the size of the simulation box
* @param inPeriodicLevel the number of level upper to 0 that will be requiried
*/
FSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inPeriodicLevel = 0)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter, inPeriodicLevel),
FSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
devM2lP(int(((inDevP*2)+1) * ((inDevP*2)+2) * 0.5)),
preM2LTransitions(0) {
allocAndInit();
......@@ -85,7 +85,7 @@ public:
/** Destructor */
~FSphericalKernel(){
if( preM2LTransitions.isLast() ){
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight + Parent::periodicLevels);
FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight);
}
}
......@@ -95,7 +95,7 @@ public:
// For all neighbors compute M2L
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
if( distantNeighbors[idxNeigh] ){
const FComplexe* const transitionVector = &preM2LTransitions[inLevel + Parent::periodicLevels][idxNeigh * devM2lP];
const FComplexe* const transitionVector = &preM2LTransitions[inLevel][idxNeigh * devM2lP];
multipoleToLocal(pole->getLocal(), distantNeighbors[idxNeigh]->getMultipole(), transitionVector);
}
}
......
Src/Kernels/Spherical/FSphericalRotationKernel.hpp
View file @ 95626445
......@@ -352,22 +352,22 @@ protected:
void allocAndInit(){
// M2L transfer, there is a maximum of 3 neighbors in each direction,
// so 6 in each dimension
preM2LTransitions = new RotationM2LTransfer*[Parent::treeHeight + Parent::periodicLevels];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe*) * (Parent::treeHeight + Parent::periodicLevels));
preM2LTransitions = new RotationM2LTransfer*[Parent::treeHeight];
memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe*) * (Parent::treeHeight));
// We start from the higher level
FReal treeWidthAtLevel = Parent::boxWidth * FReal(1 << Parent::periodicLevels);
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
FReal treeWidthAtLevel = Parent::boxWidth;
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
// Allocate data for this level
preM2LTransitions[idxLevel + Parent::periodicLevels] = reinterpret_cast<RotationM2LTransfer*>(new char[(7 * 7 * 7) * sizeof(RotationM2LTransfer)]);
preM2LTransitions[idxLevel] = reinterpret_cast<RotationM2LTransfer*>(new char[(7 * 7 * 7) * sizeof(RotationM2LTransfer)]);
// Precompute transfer vector
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
new (&preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)]) RotationM2LTransfer(Parent::devP,devM2lP,Parent::harmonic.getExpSize());
new (&preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)]) RotationM2LTransfer(Parent::devP,devM2lP,Parent::harmonic.getExpSize());
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const FPoint relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)].transfer_M2L_rotation_Fill(FSpherical(relativePos), rotation_Info);
preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)].transfer_M2L_rotation_Fill(FSpherical(relativePos), rotation_Info);
}
}
}
......@@ -385,8 +385,8 @@ public:
* @param inBoxWidth the size of the simulation box
* @param inPeriodicLevel the number of level upper to 0 that will be requiried
*/
FSphericalRotationKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const int inPeriodicLevel = 0)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter, inPeriodicLevel),
FSphericalRotationKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
devM2lP(int(((inDevP*2)+1) * ((inDevP*2)+2) * 0.5)),
preM2LTransitions(0),
rotation_Info(inDevP) {
......@@ -404,11 +404,11 @@ public:
/** Destructor */
~FSphericalRotationKernel(){
if( preM2LTransitions.isLast() ){
for(int idxLevel = -Parent::periodicLevels ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
for(int idx = 0 ; idx < 7*7*7 ; ++idx ){
preM2LTransitions[idxLevel + Parent::periodicLevels][idx].~RotationM2LTransfer();
preM2LTransitions[idxLevel][idx].~RotationM2LTransfer();
}
delete[] reinterpret_cast<char*>(preM2LTransitions[idxLevel + Parent::periodicLevels]);
delete[] reinterpret_cast<char*>(preM2LTransitions[idxLevel]);
}
}
}
......@@ -419,7 +419,7 @@ public:
// For all neighbors compute M2L
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){