Commit 642e7408 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
parents 49ecf596 bfb55710
......@@ -30,8 +30,8 @@
* Please read the license
*
* Load a file with a format like :
* NB_particles Box_width Box_X Box_Y Box_Z // init
* X Y Z // one particle by line
* NB_particles Box_width Center_X Center_Y Center_Z // init
* X Y Z PhysicalValue // one particle by line
* ....
* <code>
* FFmaBinLoader<FBasicParticle> loader("../Adir/Tests/particles.basic.txt"); <br>
......
......@@ -30,8 +30,8 @@
* Please read the license
*
* Load a file with a format like :
* NB_particles Box_width Box_X Box_Y Box_Z // init
* X Y Z // one particle by line
* NB_particles Box_width Center_X Center_Y Center_Z // init
* X Y Z PhysicalValue // one particle by line
* ....
* @code
* FFmaLoader<FBasicParticle> loader("../ADir/Tests/particles.basic.txt"); <br>
......
......@@ -45,7 +45,7 @@ public:
FRandomLoader(const size_t inNbParticles, const FReal inBoxWidth = 1.0,
const FPoint& inCenterOfBox = FPoint(0,0,0), const unsigned int inSeed = static_cast<unsigned int>(0))
: nbParticles(inNbParticles), boxWidth(inBoxWidth), centerOfBox(inCenterOfBox) {
srand(inSeed);
srand48(inSeed);
}
/**
......@@ -99,7 +99,7 @@ public:
/** Get a random number between 0 & 1 */
FReal getRandom() const{
return FReal(rand())/FReal(RAND_MAX);
return FReal(drand48());
}
};
......
......@@ -47,7 +47,8 @@ class FChebInterpolator : FNoCopyable
// compile time constants and types
enum {nnodes = TensorTraits<ORDER>::nnodes,
nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS};
nLhs = MatrixKernelClass::NLHS,
nPV = MatrixKernelClass::NPV};
typedef FChebRoots< ORDER> BasisType;
typedef FChebTensor<ORDER> TensorType;
......@@ -861,8 +862,15 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
}
}
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
FReal*const potentials = inParticles->getPotentials(idxLhs);
for(int idxLhs = 0 ; idxLhs < nLhs ; ++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
FReal targetValue = potentials[idxPart];
......@@ -1110,6 +1118,8 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
}
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPV = idxLhs % nPV;
// scale forces
forces[idxLhs][0] *= jacobian[0] / nnodes;
......@@ -1117,10 +1127,10 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
forces[idxLhs][2] *= jacobian[2] / nnodes;
// get pointers to PhysValues and force components
const FReal*const physicalValues = inParticles->getPhysicalValues(idxLhs);
FReal*const forcesX = inParticles->getForcesX(idxLhs);
FReal*const forcesY = inParticles->getForcesY(idxLhs);
FReal*const forcesZ = inParticles->getForcesZ(idxLhs);
const FReal*const physicalValues = inParticles->getPhysicalValues(idxPV);
FReal*const forcesX = inParticles->getForcesX(idxPot);
FReal*const forcesY = inParticles->getForcesY(idxPot);
FReal*const forcesZ = inParticles->getForcesZ(idxPot);
// set computed forces
forcesX[idxPart] += forces[idxLhs][0] * physicalValues[idxPart];
......@@ -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
} // 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
const FReal*const physicalValues = inParticles->getPhysicalValues(idxLhs);
FReal*const forcesX = inParticles->getForcesX(idxLhs);
FReal*const forcesY = inParticles->getForcesY(idxLhs);
FReal*const forcesZ = inParticles->getForcesZ(idxLhs);
FReal*const potentials = inParticles->getPotentials(idxLhs);
const FReal*const physicalValues = inParticles->getPhysicalValues(idxPV);
FReal*const forcesX = inParticles->getForcesX(idxPot);
FReal*const forcesY = inParticles->getForcesY(idxPot);
FReal*const forcesZ = inParticles->getForcesZ(idxPot);
FReal*const potentials = inParticles->getPotentials(idxPot);
// set computed potential
potentials[idxPart] += (potential[idxLhs]); // 1 flop
......
......@@ -56,9 +56,9 @@ class FChebSymTensorialKernel
: public FAbstractChebKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
typedef FAbstractChebKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> AbstractBaseClass;
typedef SymmetryHandler<ORDER, MatrixKernelClass::DIM, MatrixKernelClass::Type> SymmetryHandlerClass;
typedef SymmetryHandler<ORDER, MatrixKernelClass::NCMP, MatrixKernelClass::Type> SymmetryHandlerClass;
enum {nnodes = AbstractBaseClass::nnodes,
dim = MatrixKernelClass::DIM,
ncmp = MatrixKernelClass::NCMP,
nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS};
......
......@@ -512,11 +512,11 @@ static void precompute(const int idxK, const FReal CellWidth,
because it allows us to use to associate the far-field interactions based on
the index \f$t = 7^2(i+3) + 7(j+3) + (k+3)\f$ where \f$(i,j,k)\f$ denotes
the relative position of the source cell to the target cell. */
template <int ORDER, int DIM, KERNEL_FUNCTION_TYPE TYPE> class SymmetryHandler;
template <int ORDER, int NCMP, KERNEL_FUNCTION_TYPE TYPE> class SymmetryHandler;
/*! Specialization for homogeneous kernel functions */
template <int ORDER, int DIM>
class SymmetryHandler<ORDER, DIM, HOMOGENEOUS>
template <int ORDER, int NCMP>
class SymmetryHandler<ORDER, NCMP, HOMOGENEOUS>
{
static const unsigned int nnodes = ORDER*ORDER*ORDER;
......@@ -539,10 +539,10 @@ public:
{
// For each component of kernel matrix
// init all 343 item to zero, because effectively only 16 exist
K = new FReal** [DIM];
LowRank = new int* [DIM];
K = new FReal** [NCMP];
LowRank = new int* [NCMP];
for (unsigned int d=0; d<DIM; ++d) {
for (unsigned int d=0; d<NCMP; ++d) {
K[d] = new FReal* [343];
LowRank[d] = new int [343];
......@@ -565,7 +565,7 @@ public:
// precompute 16 M2L operators
const FReal ReferenceCellWidth = FReal(2.);
for (unsigned int d=0; d<DIM; ++d)
for (unsigned int d=0; d<NCMP; ++d)
precompute<ORDER,MatrixKernelClass>(d, ReferenceCellWidth, Epsilon, K[d], LowRank[d]);
}
......@@ -575,7 +575,7 @@ public:
/** Destructor */
~SymmetryHandler()
{
for (unsigned int d=0; d<DIM; ++d)
for (unsigned int d=0; d<NCMP; ++d)
for (unsigned int t=0; t<343; ++t) if (K[d][t]!=NULL) delete [] K[d][t];
}
......@@ -598,8 +598,8 @@ public:
/*! Specialization for non-homogeneous kernel functions */
template <int ORDER, int DIM>
class SymmetryHandler<ORDER, DIM, NON_HOMOGENEOUS>
template <int ORDER, int NCMP>
class SymmetryHandler<ORDER, NCMP, NON_HOMOGENEOUS>
{
static const unsigned int nnodes = ORDER*ORDER*ORDER;
......@@ -625,9 +625,9 @@ public:
: TreeHeight(inTreeHeight)
{
// init all 343 item to zero, because effectively only 16 exist
K = new FReal*** [DIM];
LowRank = new int** [DIM];
for (unsigned int d=0; d<DIM; ++d) {
K = new FReal*** [NCMP];
LowRank = new int** [NCMP];
for (unsigned int d=0; d<NCMP; ++d) {
K[d] = new FReal** [TreeHeight];
LowRank[d] = new int* [TreeHeight];
K[d][0] = NULL; K[d][1] = NULL;
......@@ -654,7 +654,7 @@ public:
}
// precompute 16 M2L operators at all levels having far-field interactions
for (unsigned int d=0; d<DIM; ++d) {
for (unsigned int d=0; d<NCMP; ++d) {
FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
CellWidth /= FReal(2.); // at level 2
for (unsigned int l=2; l<TreeHeight; ++l) {
......@@ -669,7 +669,7 @@ public:
/** Destructor */
~SymmetryHandler()
{
for (unsigned int d=0; d<DIM; ++d) {
for (unsigned int d=0; d<NCMP; ++d) {
for (unsigned int l=0; l<TreeHeight; ++l) {
if (K[d][l]!=NULL) {
for (unsigned int t=0; t<343; ++t) if (K[d][l][t]!=NULL) delete [] K[d][l][t];
......
......@@ -46,7 +46,9 @@ class FChebTensorialKernel
: public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
enum {nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS};
nLhs = MatrixKernelClass::NLHS,
nPot = MatrixKernelClass::NPOT,
nPv = MatrixKernelClass::NPV};
protected://PB: for OptiDis
......@@ -137,19 +139,19 @@ public:
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
// update local index
int idxLoc = idxV*nLhs + idxLhs;
// update local index
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
int idxMul = idxV*nRhs + idxRhs;
// update kernel index such that: x_i = K_{ij}y_j
int idxK = idxLhs*nRhs + idxRhs;
const int idxMul = idxV*nRhs + idxRhs;
// get index in matrix kernel
unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxK);
const unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
......@@ -158,8 +160,7 @@ public:
CompressedLocalExpansion);
}
}
}// NRHS
}// NLHS
}// NLHS=NPOT*NPV
}// NVALS
}
......@@ -223,14 +224,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
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*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
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
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
dim = MatrixKernelClass::DIM,
nidx = MatrixKernelClass::NIDX};
ncmp = MatrixKernelClass::NCMP};
// const MatrixKernelClass MatrixKernel;
......@@ -119,18 +118,18 @@ public:
if (U||B) throw std::runtime_error("U or B operator already set");
// allocate C
C = new FReal*[dim];
for (unsigned int d=0; d<dim; ++d)
C = new FReal*[ncmp];
for (unsigned int d=0; d<ncmp; ++d)
C[d]=NULL;
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
if (C[d]) throw std::runtime_error("Compressed M2L operator already set");
// Compute matrix of interactions
const FReal ReferenceCellWidth = FReal(2.);
rank = ComputeAndCompress<order>(MatrixKernel, ReferenceCellWidth, epsilon, U, C, B);
unsigned long sizeM2L = 343*dim*rank*rank*sizeof(FReal);
unsigned long sizeM2L = 343*ncmp*rank*rank*sizeof(FReal);
// write info
......@@ -142,7 +141,7 @@ public:
{
if (U != NULL) delete [] U;
if (B != NULL) delete [] B;
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
if (C[d] != NULL) delete [] C[d];
}
......@@ -205,8 +204,7 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCop
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
dim = MatrixKernelClass::DIM,
nidx = MatrixKernelClass::NIDX};
ncmp = MatrixKernelClass::NCMP};
// Tensorial MatrixKernel and homogeneity specific
FReal **U, **B;
......@@ -251,14 +249,14 @@ public:
// allocate C
C = new FReal**[TreeHeight];
for (unsigned int l=0; l<TreeHeight; ++l){
C[l] = new FReal*[dim];
for (unsigned int d=0; d<dim; ++d)
C[l] = new FReal*[ncmp];
for (unsigned int d=0; d<ncmp; ++d)
C[l][d]=NULL;
}
for (unsigned int l=0; l<TreeHeight; ++l) {
if (U[l] || B[l]) throw std::runtime_error("Operator U or B already set");
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
if (C[l][d]) throw std::runtime_error("Compressed M2L operator already set");
}
......@@ -270,7 +268,7 @@ public:
rank[l] = ComputeAndCompress<order>(MatrixKernel, CellWidth, epsilon, U[l], C[l], B[l]);
CellWidth /= FReal(2.); // at level l+1
}
unsigned long sizeM2L = (TreeHeight-2)*343*dim*rank[2]*rank[2]*sizeof(FReal);
unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*rank[2]*rank[2]*sizeof(FReal);
// write info
std::cout << "Compute and Set M2L operators of " << TreeHeight-2 << " levels ("<< long(sizeM2L/**1e-6*/) <<" Bytes) in "
......@@ -283,7 +281,7 @@ public:
for (unsigned int l=0; l<TreeHeight; ++l) {
if (U[l] != NULL) delete [] U[l];
if (B[l] != NULL) delete [] B[l];
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
if (C[l][d] != NULL) delete [] C[l][d];
}
}
......@@ -365,10 +363,10 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const unsigned int order = ORDER;
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
const unsigned int ninteractions = 316;
const unsigned int dim = MatrixKernelClass::DIM;
const unsigned int ncmp = MatrixKernelClass::NCMP;
// allocate memory and store compressed M2L operators
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
if (C[d]) throw std::runtime_error("Compressed M2L operators are already set");
// interpolation points of source (Y) and target (X) cell
......@@ -378,8 +376,8 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
// allocate memory and compute 316 m2l operators
FReal** _C;
_C = new FReal* [dim];
for (unsigned int d=0; d<dim; ++d)
_C = new FReal* [ncmp];
for (unsigned int d=0; d<ncmp; ++d)
_C[d] = new FReal [nnodes*nnodes * ninteractions];
unsigned int counter = 0;
......@@ -398,11 +396,11 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
// Compute current M2L interaction (block matrix)
FReal* block;
block = new FReal[dim];
block = new FReal[ncmp];
MatrixKernel->evaluateBlock(X[m], Y[n], block);
// Copy block in C
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
_C[d][counter*nnodes*nnodes + n*nnodes + m] = block[d];
delete [] block;
......@@ -446,24 +444,24 @@ unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
// store C
counter = 0;
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
C[d] = new FReal [343 * rank*rank];
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
FBlas::copy(rank*rank, _C[d] + counter*rank*rank, C[d] + idx*rank*rank);
counter++;
} else {
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
FBlas::setzero(rank*rank, C[d] + idx*rank*rank);
}
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
for (unsigned int d=0; d<dim; ++d)
for (unsigned int d=0; d<ncmp; ++d)
delete [] _C[d];
......
......@@ -44,8 +44,8 @@ enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
* This class provides the evaluators and scaling functions of the matrix
* kernels. A matrix kernel should be understood in the sense of a kernel
* of interaction (or the fundamental solution of a given equation).
* It can either be scalar (DIM=1) or tensorial (DIM>1) depending on the
* dimension of the equation considered. DIM denotes the number of components
* It can either be scalar (NCMP=1) or tensorial (NCMP>1) depending on the
* dimension of the equation considered. NCMP denotes the number of components
* that are actually stored (e.g. 6 for a \f$3\times3\f$ symmetric tensor).
* Notes on application scheme:
* Let there be a kernel \f$K\f$ such that \f$X_i=K_{ij}Y_j\f$
......@@ -70,10 +70,11 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
static const unsigned int NCMP = 1; //< number of components
static const unsigned int NPV = 1; //< dim of physical values
static const unsigned int NPOT = 1; //< dim of potentials
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelR(const double = 0.0, const unsigned int = 0) {}
......@@ -126,10 +127,11 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
static const unsigned int NCMP = 1; //< number of components
static const unsigned int NPV = 1; //< dim of physical values
static const unsigned int NPOT = 1; //< dim of potentials
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelRR(const double = 0.0, const unsigned int = 0) {}
......@@ -145,6 +147,11 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -164,10 +171,11 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = LENNARD_JONES_POTENTIAL;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
static const unsigned int NCMP = 1; //< number of components
static const unsigned int NPV = 1; //< dim of physical values
static const unsigned int NPOT = 1; //< dim of potentials
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelLJ(const double = 0.0, const unsigned int = 0) {}
......@@ -186,6 +194,11 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
return one_over_r6 * one_over_r6 - one_over_r6;
}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
FReal getScaleFactor(const FReal, const int) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
......@@ -203,7 +216,27 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//
// Tensorial Matrix Kernels (DIM>1)
// Tensorial Matrix Kernels (NCMP>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 +246,15 @@ 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 DIM = 6; //PB: dimension of kernel
static const unsigned int NIDX = 2; //PB: number of indices
static const unsigned int indexTab[/*DIM*NIDX=12*/];
static const unsigned int NRHS = 3;
static const unsigned int NLHS = 3;
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*/];
......@@ -225,7 +262,7 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
const unsigned int _i,_j;
FInterpMatrixKernel_IOR(const double = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM])
: _i(indexTab[d]), _j(indexTab[d+NCMP])
{}
// returns position in reduced storage from position in full 3x3 matrix
......@@ -249,9 +286,9 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm();
for(unsigned int d=0;d<DIM;++d){
for(unsigned int d=0;d<NCMP;++d){
// unsigned int i = indexTab[d];
// unsigned int j = indexTab[d+DIM];
// unsigned int j = indexTab[d+NCMP];
// if(i==j)