Skip to content
GitLab
Commit 61c7edb8 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files

Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

parents 34f3ff0e 088fbae5
Hide whitespace changes
Inline Side-by-side
Src/Kernels/Interpolation/FInterpMatrixKernel.cpp
View file @ 61c7edb8
#include "FInterpMatrixKernel.hpp"
/// ID_OVER_R
const unsigned int FInterpMatrixKernel_IOR::indexTab[]={0,0,0,1,1,2,
0,1,2,1,2,2};
const unsigned int FInterpMatrixKernel_IOR::applyTab[]={0,1,2,
1,3,4,
2,4,5};
/// R_IJ
const unsigned int FInterpMatrixKernel_R_IJ::indexTab[]={0,0,0,1,1,2,
0,1,2,1,2,2};
......
Src/Kernels/Interpolation/FInterpMatrixKernel.hpp
View file @ 61c7edb8
......@@ -27,7 +27,6 @@
enum KERNEL_FUNCTION_IDENTIFIER {ONE_OVER_R,
ONE_OVER_R_SQUARED,
LENNARD_JONES_POTENTIAL,
ID_OVER_R,
R_IJ,
R_IJK};
......@@ -68,7 +67,7 @@ struct FInterpAbstractMatrixKernel : FNoCopyable
/// One over r
struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R;
static const unsigned int NCMP = 1; //< number of components
static const unsigned int NPV = 1; //< dim of physical values
......@@ -106,18 +105,6 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
return FReal(2.) / CellWidth;
}
// FReal getScaleFactor(const FReal, const int) const
// {
// // return 1 because non homogeneous kernel functions cannot be scaled!!!
// return FReal(1.);
// }
//
// FReal getScaleFactor(const FReal) const
// {
// // return 1 because non homogeneous kernel functions cannot be scaled!!!
// return FReal(1.);
// }
};
......@@ -160,7 +147,7 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(4.) / CellWidth;
return FReal(4.) / (CellWidth*CellWidth);
}
};
......@@ -241,80 +228,6 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/// Test Tensorial kernel 1/R*Id_3
struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R;
static const unsigned int NK = 3*3; //< total number of components
static const unsigned int NCMP = 6; //< number of components after symmetry
static const unsigned int NPV = 3; //< dim of physical values
static const unsigned int NPOT = 3; //< dim of potentials
static const unsigned int NRHS = NPV; //< dim of mult exp
static const unsigned int NLHS = NPOT*NRHS; //< dim of loc exp
// store indices (i,j) corresponding to sym idx
static const unsigned int indexTab[/*2*NCMP=12*/];
// store positions in sym tensor
static const unsigned int applyTab[/*9*/];
const unsigned int _i,_j;
FInterpMatrixKernel_IOR(const double = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+NCMP])
{}
// returns position in reduced storage from position in full 3x3 matrix
unsigned int getPosition(const unsigned int n) const
{return applyTab[n];}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
// low rank approx does not support nul kernels
// if(_i==_j)
return FReal(1.)/xy.norm();
// else
// return FReal(0.);
}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm();
for(unsigned int d=0;d<NCMP;++d){
// unsigned int i = indexTab[d];
// unsigned int j = indexTab[d+NCMP];
// if(i==j)
block[d] = one_over_r;
// else
// block[d] = 0.0;
}
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(2.) / CellWidth;
}
// FReal getScaleFactor(const FReal) const
// {
// // return 1 because non homogeneous kernel functions cannot be scaled!!!
// return FReal(1.);
// }
};
/// R_{,ij}
// PB: IMPORTANT! This matrix kernel does not present the symmetries
......@@ -413,12 +326,6 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
return FReal(2.) / CellWidth;
}
// // R_{,ij} is set non-homogeneous
// FReal getScaleFactor(const FReal CellWidth) const
// {
// return FReal(1.);
// }
};
/// R_{,ijk}
......@@ -550,12 +457,6 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
return FReal(4.) / (CellWidth*CellWidth);
}
// // R_{,ijk} is set non-homogeneous
// FReal getScaleFactor(const FReal CellWidth) const
// {
// return FReal(1.);
// }
};
......
Src/Kernels/Interpolation/FInterpP2PKernels.hpp
View file @ 61c7edb8
......@@ -48,24 +48,6 @@ struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, 1>
}
};
/*! Specialization for ID_OVER_R potential */
template <>
struct DirectInteractionComputer<ID_OVER_R, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualIOR(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteIOR(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for GradGradR potential */
template <>
struct DirectInteractionComputer<R_IJ, 1>
......@@ -156,27 +138,6 @@ struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, NVALS>
}
};
/*! Specialization for ID_OVER_R potential */
template <int NVALS>
struct DirectInteractionComputer<ID_OVER_R, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualIOR(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemoteIOR(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for GradGradR potential */
template <int NVALS>
......
Src/Kernels/P2P/FP2P.hpp
View file @ 61c7edb8
......@@ -1051,301 +1051,6 @@ namespace FP2P {
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// PB: Test Tensorial Kernel ID_over_R
// i.e. [[1 1 1]
// [1 1 1] * 1/R
// [1 1 1]]
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
/**
* @brief FullMutualIOR
*/
template <class ContainerClass>
inline void FullMutualIOR(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors){
const int nbParticlesTargets = inTargets->getNbParticles();
const FReal*const targetsX = inTargets->getPositions()[0];
const FReal*const targetsY = inTargets->getPositions()[1];
const FReal*const targetsZ = inTargets->getPositions()[2];
for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials(i);
FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX(i);
FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY(i);
FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ(i);
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
//if(i==j)
potentialCoef = inv_distance;
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
//}
//else{
// coef[k] = FReal(0.);
//}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
sourcesForcesX[idxSource] -= coef[0];
sourcesForcesY[idxSource] -= coef[1];
sourcesForcesZ[idxSource] -= coef[2];
sourcesPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
}
}
}
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
FReal dx = targetsX[idxSource] - targetsX[idxTarget];
FReal dy = targetsY[idxSource] - targetsY[idxTarget];
FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
// potentials
FReal potentialCoef;
//if(i==j)
potentialCoef = inv_distance;
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource]);
//}
//else{
// coef[k] = FReal(0.);
//}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * targetsPhysicalValues[idxSource] );
targetsForcesX[idxSource] -= coef[0];
targetsForcesY[idxSource] -= coef[1];
targetsForcesZ[idxSource] -= coef[2];
targetsPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
}
}
/**
* @brief FullRemoteIOR
*/
template <class ContainerClass>
inline void FullRemoteIOR(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors){
const int nbParticlesTargets = inTargets->getNbParticles();
const FReal*const targetsX = inTargets->getPositions()[0];
const FReal*const targetsY = inTargets->getPositions()[1];
const FReal*const targetsZ = inTargets->getPositions()[2];
for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
// potential
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
//if(i==j)
potentialCoef = inv_distance ;
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
//}
//else{
// coef[k] = FReal(0.);
//}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
}// j
}// i
}
}
}
}
}
/**
* @brief MutualParticlesIOR
* @param sourceX
* @param sourceY
* @param sourceZ
* @param sourcePhysicalValue
* @param sourceForceX
* @param sourceForceY
* @param sourceForceZ
* @param sourcePotential
* @param targetX
* @param targetY
* @param targetZ
* @param targetPhysicalValue
* @param targetForceX
* @param targetForceY
* @param targetForceZ
* @param targetPotential
*/
inline void MutualParticlesIOR(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal* sourcePhysicalValue,
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal* targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential
){
// GradGradR potential
FReal dx = sourceX - targetX;
FReal dy = sourceY - targetY;
FReal dz = sourceZ - targetZ;
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
for(unsigned int j = 0 ; j < 3 ; ++j){
// potentials
FReal potentialCoef;
//if(i==j)
potentialCoef = inv_distance;
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetPhysicalValue[j] * sourcePhysicalValue[j]);
//}
//else{
// coef[k] = FReal(0.);
//}
}// k
targetForceX[i] += coef[0];
targetForceY[i] += coef[1];
targetForceZ[i] += coef[2];
targetPotential[i] += ( potentialCoef * sourcePhysicalValue[j] );
sourceForceX[i] -= coef[0];
sourceForceY[i] -= coef[1];
sourceForceZ[i] -= coef[2];
sourcePotential[i] += ( potentialCoef * targetPhysicalValue[j] );
}// j
}// i
}
}
#ifdef ScalFMM_USE_SSE
......
Src/Kernels/Uniform/FUnifKernel.hpp
View file @ 61c7edb8
......@@ -52,7 +52,8 @@ class FUnifKernel
AbstractBaseClass;
/// Needed for M2L operator
FSmartPointer< M2LHandlerClass,FSmartPointerMemory> M2LHandler;
const M2LHandlerClass M2LHandler;
public:
/**
......@@ -64,9 +65,9 @@ public:
const FReal inBoxWidth,
const FPoint& inBoxCenter)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight,
inBoxWidth))// PB: for non homogeneous case
M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight,
inBoxWidth)
{ }
......@@ -79,7 +80,7 @@ public:
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxRhs), SourceParticles);
// 2) apply Discrete Fourier Transform
M2LHandler->applyZeroPaddingAndDFT(LeafCell->getMultipole(idxRhs),
M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxRhs),
LeafCell->getTransformedMultipole(idxRhs));
}
......@@ -100,7 +101,7 @@ public:
}
}
// 2) Apply Discete Fourier Transform
M2LHandler->applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs),
M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxRhs),
ParentCell->getTransformedMultipole(idxRhs));
}
}
......@@ -137,7 +138,7 @@ public:
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
M2LHandler->applyFC(idx, TreeLevel, scale,
M2LHandler.applyFC(idx, TreeLevel, scale,
SourceCells[idx]->getTransformedMultipole(idxRhs),
TransformedLocalExpansion);
}
......@@ -168,7 +169,7 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) Apply Inverse Discete Fourier Transform
M2LHandler->unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxRhs),
M2LHandler.unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxRhs),
const_cast<CellClass*>(ParentCell)->getLocal(idxRhs));
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
......@@ -187,7 +188,7 @@ public:
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) Apply Inverse Discete Fourier Transform
M2LHandler->unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs),
M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs),
const_cast<CellClass*>(LeafCell)->getLocal(idxRhs));
// 2.a) apply Sx
......
Src/Kernels/Uniform/FUnifM2LHandler.hpp
View file @ 61c7edb8
......@@ -40,7 +40,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
// allocate memory and store compressed M2L operators
if (FC) throw std::runtime_error("M2L operators are already set");
// PB: need to redefine some constant since not function from m2lhandler class
// dimensions of operators
const unsigned int order = ORDER;
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
const unsigned int ninteractions = 316;
......@@ -68,8 +68,8 @@ static void Compute(const MatrixKernelClass *const MatrixKernel, const FReal Cel
// init Discrete Fourier Transformator
const int dimfft = 1; // unidim FFT since fully circulant embedding
const int steps[dimfft] = {rc};