Commit 1c8ab4a7 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Introduced core width in R_IJ and R_IJK kernels; Fixed R_IJK (scheme:...

Introduced core width in R_IJ and R_IJK kernels; Fixed R_IJK (scheme: p_i=R_ijk*q_jk, BEWARE: antisymmetry and non_homogeneity strongly affect P2P and MatrixKernel, force not yet implemented); Improved consistency with potential theory for tensorial matrix kernels by defining number of phys vals NPV and number of potential vals NPOT (NRHS=NPV, NLHS=NPOT*NRHS).
parent 2f6e90c4
...@@ -47,7 +47,8 @@ class FChebInterpolator : FNoCopyable ...@@ -47,7 +47,8 @@ class FChebInterpolator : FNoCopyable
// compile time constants and types // compile time constants and types
enum {nnodes = TensorTraits<ORDER>::nnodes, enum {nnodes = TensorTraits<ORDER>::nnodes,
nRhs = MatrixKernelClass::NRHS, nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS}; nLhs = MatrixKernelClass::NLHS,
nPV = MatrixKernelClass::NPV};
typedef FChebRoots< ORDER> BasisType; typedef FChebRoots< ORDER> BasisType;
typedef FChebTensor<ORDER> TensorType; typedef FChebTensor<ORDER> TensorType;
...@@ -861,8 +862,15 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c ...@@ -861,8 +862,15 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
} }
} }
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
FReal*const potentials = inParticles->getPotentials(idxLhs); // distribution over potential components:
// We sum the multidim contribution of PhysValue
// This was originally done at M2L step but moved here
// because their storage is required by the force computation.
// In fact : f_{ik}(x)=w_j(x) \nabla_{x_i} K_{ij}(x,y)w_j(y))
const unsigned int idxPot = idxLhs / nPV;
FReal*const potentials = inParticles->getPotentials(idxPot);
// interpolate and increment target value // interpolate and increment target value
FReal targetValue = potentials[idxPart]; FReal targetValue = potentials[idxPart];
...@@ -1110,6 +1118,8 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -1110,6 +1118,8 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
} }
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPV = idxLhs % nPV;
// scale forces // scale forces
forces[idxLhs][0] *= jacobian[0] / nnodes; forces[idxLhs][0] *= jacobian[0] / nnodes;
...@@ -1117,10 +1127,10 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -1117,10 +1127,10 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
forces[idxLhs][2] *= jacobian[2] / nnodes; forces[idxLhs][2] *= jacobian[2] / nnodes;
// get pointers to PhysValues and force components // get pointers to PhysValues and force components
const FReal*const physicalValues = inParticles->getPhysicalValues(idxLhs); const FReal*const physicalValues = inParticles->getPhysicalValues(idxPV);
FReal*const forcesX = inParticles->getForcesX(idxLhs); FReal*const forcesX = inParticles->getForcesX(idxPot);
FReal*const forcesY = inParticles->getForcesY(idxLhs); FReal*const forcesY = inParticles->getForcesY(idxPot);
FReal*const forcesZ = inParticles->getForcesZ(idxLhs); FReal*const forcesZ = inParticles->getForcesZ(idxPot);
// set computed forces // set computed forces
forcesX[idxPart] += forces[idxLhs][0] * physicalValues[idxPart]; forcesX[idxPart] += forces[idxLhs][0] * physicalValues[idxPart];
...@@ -1342,12 +1352,15 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi ...@@ -1342,12 +1352,15 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
forces[idxLhs][2] = ( FReal(2.)*f2[3] + FReal(4.)*f4[3] + FReal(8.)*f8[3]) * jacobian[2] / nnodes; // 7 flops forces[idxLhs][2] = ( FReal(2.)*f2[3] + FReal(4.)*f4[3] + FReal(8.)*f8[3]) * jacobian[2] / nnodes; // 7 flops
} // 28 + (ORDER-1) * ((ORDER-1) * (27 + (ORDER-1) * 16)) flops } // 28 + (ORDER-1) * ((ORDER-1) * (27 + (ORDER-1) * 16)) flops
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPV = idxLhs % nPV;
// get potentials, physValues and forces components // get potentials, physValues and forces components
const FReal*const physicalValues = inParticles->getPhysicalValues(idxLhs); const FReal*const physicalValues = inParticles->getPhysicalValues(idxPV);
FReal*const forcesX = inParticles->getForcesX(idxLhs); FReal*const forcesX = inParticles->getForcesX(idxPot);
FReal*const forcesY = inParticles->getForcesY(idxLhs); FReal*const forcesY = inParticles->getForcesY(idxPot);
FReal*const forcesZ = inParticles->getForcesZ(idxLhs); FReal*const forcesZ = inParticles->getForcesZ(idxPot);
FReal*const potentials = inParticles->getPotentials(idxLhs); FReal*const potentials = inParticles->getPotentials(idxPot);
// set computed potential // set computed potential
potentials[idxPart] += (potential[idxLhs]); // 1 flop potentials[idxPart] += (potential[idxLhs]); // 1 flop
......
...@@ -46,7 +46,9 @@ class FChebTensorialKernel ...@@ -46,7 +46,9 @@ class FChebTensorialKernel
: public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> : public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{ {
enum {nRhs = MatrixKernelClass::NRHS, enum {nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS}; nLhs = MatrixKernelClass::NLHS,
nPot = MatrixKernelClass::NPOT,
nPv = MatrixKernelClass::NPV};
protected://PB: for OptiDis protected://PB: for OptiDis
...@@ -137,19 +139,19 @@ public: ...@@ -137,19 +139,19 @@ public:
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){ for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
// update local index // update local index
int idxLoc = idxV*nLhs + idxLhs; const int idxLoc = idxV*nLhs + idxLhs;
FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxLoc) + AbstractBaseClass::nnodes; FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxLoc) + AbstractBaseClass::nnodes;
for (int idxRhs=0; idxRhs < nRhs; ++idxRhs){ // update idxRhs
const int idxRhs = idxLhs % nPv;
// update multipole index // update multipole index
int idxMul = idxV*nRhs + idxRhs; const int idxMul = idxV*nRhs + idxRhs;
// update kernel index such that: x_i = K_{ij}y_j
int idxK = idxLhs*nRhs + idxRhs;
// get index in matrix kernel // get index in matrix kernel
unsigned int d const unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxK); = AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
for (int idx=0; idx<343; ++idx){ for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){ if (SourceCells[idx]){
...@@ -158,8 +160,7 @@ public: ...@@ -158,8 +160,7 @@ public:
CompressedLocalExpansion); CompressedLocalExpansion);
} }
} }
}// NRHS }// NLHS=NPOT*NPV
}// NLHS
}// NVALS }// NVALS
} }
...@@ -223,14 +224,14 @@ public: ...@@ -223,14 +224,14 @@ public:
ContainerClass* const NeighborSourceParticles[27], ContainerClass* const NeighborSourceParticles[27],
const int /* size */) const int /* size */)
{ {
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles); DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
} }
void P2PRemote(const FTreeCoordinate& /*inPosition*/, void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/, ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/){ ContainerClass* const inNeighbors[27], const int /*inSize*/){
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27); DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
} }
}; };
......
...@@ -86,8 +86,7 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl ...@@ -86,8 +86,7 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
enum {order = ORDER, enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes, nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field) ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
dim = MatrixKernelClass::DIM, dim = MatrixKernelClass::DIM};
nidx = MatrixKernelClass::NIDX};
// const MatrixKernelClass MatrixKernel; // const MatrixKernelClass MatrixKernel;
...@@ -205,8 +204,7 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop ...@@ -205,8 +204,7 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
enum {order = ORDER, enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes, nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field) ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
dim = MatrixKernelClass::DIM, dim = MatrixKernelClass::DIM};
nidx = MatrixKernelClass::NIDX};
// Tensorial MatrixKernel and homogeneity specific // Tensorial MatrixKernel and homogeneity specific
FReal **U, **B; FReal **U, **B;
......
...@@ -71,7 +71,8 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel ...@@ -71,7 +71,8 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R; static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R;
static const unsigned int DIM = 1; //PB: dimension of kernel static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices static const unsigned int NPV = 1;
static const unsigned int NPOT = 1;
static const unsigned int NRHS = 1; static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1; static const unsigned int NLHS = 1;
...@@ -127,7 +128,8 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel ...@@ -127,7 +128,8 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED; static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
static const unsigned int DIM = 1; //PB: dimension of kernel static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices static const unsigned int NPV = 1;
static const unsigned int NPOT = 1;
static const unsigned int NRHS = 1; static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1; static const unsigned int NLHS = 1;
...@@ -165,7 +167,8 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel ...@@ -165,7 +167,8 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = LENNARD_JONES_POTENTIAL; static const KERNEL_FUNCTION_IDENTIFIER Identifier = LENNARD_JONES_POTENTIAL;
static const unsigned int DIM = 1; //PB: dimension of kernel static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices static const unsigned int NPV = 1;
static const unsigned int NPOT = 1;
static const unsigned int NRHS = 1; static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1; static const unsigned int NLHS = 1;
...@@ -205,6 +208,26 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel ...@@ -205,6 +208,26 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
// //
// Tensorial Matrix Kernels (DIM>1) // Tensorial Matrix Kernels (DIM>1)
// //
// The definition of the potential p and force f are extended to the case
// of tensorial interaction kernels:
// p_i(x) = K_{ip}(x,y)w_p(y), \forall i=1..NPOT, p=1..NPV
// f_{ik}= w_p(x)K_{ip,k}(x,y)w_p(y) "
//
// Since the interpolation scheme is such that
// p_i(x) \approx S^m(x) L^{m}_{ip}
// f_{ik}= w_p(x) \nabla_k S^m(x) L^{m}_{ip}
// with
// L^{m}_{ip} = K^{mn}_{ip} S^n(y) w_p(y) (local expansion)
// M^{m}_{p} = S^n(y) w_p(y) (multipole expansion)
// then the multipole exp have NPV components and the local exp NPOT*NPV.
//
// NB1: Only the computation of forces requires that the sum over p is
// performed at L2P step. It could be done at M2L step for the potential.
//
// NB2: An efficient application of the matrix kernel is highly kernel
// dependent, we recommand overriding the P2M/M2L/L2P function of the kernel
// you are using in order to have optimal performances + set your own NRHS/NLHS.
//
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
...@@ -213,11 +236,15 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel ...@@ -213,11 +236,15 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
{ {
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R; static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R;
static const unsigned int DIM = 6; //PB: dimension of kernel static const unsigned int NK = 3*3; // dim of kernel
static const unsigned int NIDX = 2; //PB: number of indices static const unsigned int DIM = 6; // reduced dim
static const unsigned int indexTab[/*DIM*NIDX=12*/]; static const unsigned int NPOT = 3; // dim of potentials
static const unsigned int NRHS = 3; static const unsigned int NPV = 3; // dim of physical values
static const unsigned int NLHS = 3; static const unsigned int NRHS = NPV; // dim of mult exp (here = dim physval)
static const unsigned int NLHS = NPOT*NRHS; // dim of loc exp (here = NPOT*NRHS)
// store indices (i,j) corresponding to sym idx
static const unsigned int indexTab[/*2*DIM=12*/];
// store positions in sym tensor // store positions in sym tensor
static const unsigned int applyTab[/*9*/]; static const unsigned int applyTab[/*9*/];
...@@ -271,12 +298,6 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel ...@@ -271,12 +298,6 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
return FReal(2.) / CellWidth; 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 // FReal getScaleFactor(const FReal) const
// { // {
// // return 1 because non homogeneous kernel functions cannot be scaled!!! // // return 1 because non homogeneous kernel functions cannot be scaled!!!
...@@ -290,31 +311,45 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel ...@@ -290,31 +311,45 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
// required by ChebSym kernel => only suited for Unif interpolation // required by ChebSym kernel => only suited for Unif interpolation
struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
{ {
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJ; static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJ;
static const unsigned int DIM = 6; //PB: dimension of kernel static const unsigned int NK = 3*3; // dim of kernel
static const unsigned int NIDX = 2; //PB: number of indices static const unsigned int DIM = 6; // reduced dim
static const unsigned int indexTab[/*DIM*NIDX=12*/]; static const unsigned int NPOT = 3; // dim of potentials
static const unsigned int NRHS = 3; static const unsigned int NPV = 3; // dim of physical values
static const unsigned int NLHS = 3; static const unsigned int NRHS = NPV; // dim of mult exp (here = dim physval)
static const unsigned int NLHS = NPOT*NRHS; // dim of loc exp (here = NPOT*NRHS)
// store indices (i,j) corresponding to sym idx
static const unsigned int indexTab[/*2*DIM=12*/];
// store positions in sym tensor (when looping over NRHSxNLHS) // store positions in sym tensor (when looping over NRHSxNLHS)
static const unsigned int applyTab[/*9*/]; static const unsigned int applyTab[/*NK=9*/];
// indices to be set at construction if component-wise evaluation is performed
const unsigned int _i,_j; const unsigned int _i,_j;
FInterpMatrixKernel_R_IJ(const double = 0.0, const unsigned int d = 0) // Material Parameters
: _i(indexTab[d]), _j(indexTab[d+DIM]) const double _CoreWidth2; // if >0 then kernel is NON homogeneous
FInterpMatrixKernel_R_IJ(const double CoreWidth = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM]), _CoreWidth2(CoreWidth*CoreWidth)
{} {}
// returns position in reduced storage from position in full 3x3 matrix // returns position in reduced storage from position in full 3x3 matrix
int getPosition(const unsigned int n) const int getPosition(const unsigned int n) const
{return applyTab[n];} {return applyTab[n];}
// returns Core Width squared
double getCoreWidth2() const
{return _CoreWidth2;}
FReal evaluate(const FPoint& x, const FPoint& y) const FReal evaluate(const FPoint& x, const FPoint& y) const
{ {
const FPoint xy(x-y); const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm(); const FReal one_over_r = FReal(1.)/sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ() + _CoreWidth2);
const FReal one_over_r3 = one_over_r*one_over_r*one_over_r; const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
double ri,rj; double ri,rj;
...@@ -338,7 +373,9 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel ...@@ -338,7 +373,9 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{ {
const FPoint xy(x-y); const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm(); const FReal one_over_r = FReal(1.)/sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ() + _CoreWidth2);
const FReal one_over_r3 = one_over_r*one_over_r*one_over_r; const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
const double r[3] = {xy.getX(),xy.getY(),xy.getZ()}; const double r[3] = {xy.getX(),xy.getY(),xy.getZ()};
...@@ -366,36 +403,57 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel ...@@ -366,36 +403,57 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
return FReal(2.) / CellWidth; return FReal(2.) / CellWidth;
} }
// // R_{,ij} is set non-homogeneous
// FReal getScaleFactor(const FReal CellWidth) const
// {
// return FReal(1.);
// }
}; };
/// R_{,ijk} /// R_{,ijk}
struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
{ {
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS; static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJK; static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJK;
static const unsigned int DIM = 10; //PB: dimension of kernel static const unsigned int NK = 3*3*3; // dim of kernel
static const unsigned int NIDX = 3; //PB: number of indices static const unsigned int DIM = 10; // reduced dim
static const unsigned int indexTab[/*DIM*NIDX=30*/]; static const unsigned int NPOT = 3; // dim of potentials
static const unsigned int NRHS = 3; static const unsigned int NPV = 3*3; // dim of physical values
static const unsigned int NLHS = 3; static const unsigned int NRHS = NPV; // dim of mult exp (here = dim physval)
static const unsigned int NLHS = NPOT*NRHS; // dim of loc exp (here = NPOT*NRHS)
// store positions in sym tensor (when looping over NRHSxNLHS) // store indices (i,j,k) corresponding to sym idx
static const unsigned int applyTab[/*27*/]; static const unsigned int indexTab[/*3*DIM=30*/];
// store positions in sym tensor wr to full
static const unsigned int applyTab[/*NK=27*/];
// indices to be set at construction if component-wise evaluation is performed
const unsigned int _i,_j,_k; const unsigned int _i,_j,_k;
FInterpMatrixKernel_R_IJK(const double = 0.0, const unsigned int d = 0) // Material Parameters
: _i(indexTab[d]), _j(indexTab[d+DIM]), _k(indexTab[d+2*DIM]) const double _CoreWidth2; // if >0 then kernel is NON homogeneous
FInterpMatrixKernel_R_IJK(const double CoreWidth = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM]), _k(indexTab[d+2*DIM]), _CoreWidth2(CoreWidth*CoreWidth)
{} {}
// returns position in reduced storage from position in full 3x3x3 matrix // returns position in reduced storage from position in full 3x3x3 matrix
int getPosition(const unsigned int n) const int getPosition(const unsigned int n) const
{return applyTab[n];} {return applyTab[n];}
// returns Core Width squared
double getCoreWidth2() const
{return _CoreWidth2;}
FReal evaluate(const FPoint& x, const FPoint& y) const FReal evaluate(const FPoint& x, const FPoint& y) const
{ {
const FPoint xy(x-y); // const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm(); const FPoint xy(y-x);
const FReal one_over_r = FReal(1.)/sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ() + _CoreWidth2);
const FReal one_over_r2 = one_over_r*one_over_r; const FReal one_over_r2 = one_over_r*one_over_r;
const FReal one_over_r3 = one_over_r2*one_over_r; const FReal one_over_r3 = one_over_r2*one_over_r;
double ri,rj,rk; double ri,rj,rk;
...@@ -434,6 +492,39 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel ...@@ -434,6 +492,39 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
} }
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
// const FPoint xy(x-y);
const FPoint xy(y-x);
const FReal one_over_r = FReal(1.)/sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ() + _CoreWidth2);
const FReal one_over_r2 = one_over_r*one_over_r;
const FReal one_over_r3 = one_over_r2*one_over_r;
const double r[3] = {xy.getX(),xy.getY(),xy.getZ()};
for(unsigned int d=0;d<DIM;++d){
unsigned int i = indexTab[d];
unsigned int j = indexTab[d+DIM];
unsigned int k = indexTab[d+2*DIM];
if(i==j){
if(j==k) //i=j=k
block[d] = FReal(3.) * ( FReal(-1.) + r[i]*r[i] * one_over_r2 ) * r[i] * one_over_r3;
else //i=j!=k
block[d] = ( FReal(-1.) + FReal(3.) * r[i]*r[i] * one_over_r2 ) * r[k] * one_over_r3;
}
else //(i!=j)
if(i==k) //i=k!=j
block[d] = ( FReal(-1.) + FReal(3.) * r[i]*r[i] * one_over_r2 ) * r[j] * one_over_r3;
else if(j==k) //i!=k=j
block[d] = ( FReal(-1.) + FReal(3.) * r[j]*r[j] * one_over_r2 ) * r[i] * one_over_r3;
else //i!=k!=j
block[d] = FReal(3.) * r[i]*r[j]*r[k] * one_over_r2 * one_over_r3;
}
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{ {
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel))); const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
...@@ -446,6 +537,13 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel ...@@ -446,6 +537,13 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
{ {
return FReal(4.) / (CellWidth*CellWidth); return FReal(4.) / (CellWidth*CellWidth);
} }
// // R_{,ijk} is set non-homogeneous
// FReal getScaleFactor(const FReal CellWidth) const
// {
// return FReal(1.);
// }
}; };
......
...@@ -70,17 +70,39 @@ struct DirectInteractionComputer<ID_OVER_R, 1> ...@@ -70,17 +70,39 @@ struct DirectInteractionComputer<ID_OVER_R, 1>
template <> template <>
struct DirectInteractionComputer<R_IJ, 1> struct DirectInteractionComputer<R_IJ, 1>
{ {
template <typename ContainerClass> template <typename ContainerClass, typename MatrixKernelClass>
static void P2P( ContainerClass* const FRestrict TargetParticles, static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){ ContainerClass* const NeighborSourceParticles[27],
FP2P::FullMutualRIJ(TargetParticles,NeighborSourceParticles,14);