Commit d046dd55 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Modified interpolation based kernels such that matrix kernel is passed to ctor.

parent 72b93562
......@@ -122,6 +122,7 @@ int main(int argc, char* argv[])
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//
#ifdef _OPENMP
......@@ -161,7 +162,7 @@ int main(int argc, char* argv[])
time.tic();
//
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
//
FmmClass algorithm(&tree, &kernels);
//
......
......@@ -49,8 +49,6 @@ protected:
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for P2P operator
/*const*/ FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Height of the entire oct-tree
const unsigned int TreeHeight;
/// Corner of oct-tree box
......@@ -61,9 +59,6 @@ protected:
const FReal BoxWidthLeaf;
/// Extension of the box width ( same for all level! )
const FReal BoxWidthExtension;
/// Parameter to pass to matrix kernel (material specific or anything)
const FReal MatParam;
/**
* Compute center of leaf cell from its tree coordinate.
......@@ -86,18 +81,15 @@ public:
FAbstractChebKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal inBoxWidthExtension = 0.0,
const FReal inMatParam = 0.0)
const FReal inBoxWidthExtension = 0.0)
: Interpolator(new InterpolatorClass(inTreeHeight,
inBoxWidth,
inBoxWidthExtension)),
MatrixKernel(new MatrixKernelClass(inMatParam)),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
BoxWidthExtension(inBoxWidthExtension),
MatParam(inMatParam)
BoxWidthExtension(inBoxWidthExtension)
{
/* empty */
}
......@@ -109,10 +101,6 @@ public:
const InterpolatorClass * getPtrToInterpolator() const
{ return Interpolator.getPtr(); }
MatrixKernelClass * getPtrMatrixKernel()
{return MatrixKernel; }
virtual void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles) = 0;
......
......@@ -51,6 +51,9 @@ class FChebKernel
typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for M2L operator
FSmartPointer< M2LHandlerClass,FSmartPointerMemory> M2LHandler;
......@@ -60,20 +63,21 @@ public:
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter,
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxWidth,
inBoxCenter),
M2LHandler(new M2LHandlerClass(Epsilon))
FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel,
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxWidth,
inBoxCenter),
MatrixKernel(inMatrixKernel),
M2LHandler(new M2LHandlerClass(MatrixKernel,Epsilon))
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet();
M2LHandler->ComputeAndCompressAndSet();
}
FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
: FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
: FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
{
}
......@@ -212,14 +216,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
}
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,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
}
};
......
......@@ -70,7 +70,7 @@ class FChebM2LHandler : FNoCopyable
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)
const MatrixKernelClass MatrixKernel;
const MatrixKernelClass *const MatrixKernel;
FReal *U, *C, *B;
const FReal epsilon; //<! accuracy which determines trucation of SVD
......@@ -88,8 +88,8 @@ class FChebM2LHandler : FNoCopyable
public:
FChebM2LHandler(const FReal _epsilon)
: MatrixKernel(), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
FChebM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal _epsilon)
: MatrixKernel(inMatrixKernel), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
{}
~FChebM2LHandler()
......@@ -108,7 +108,7 @@ public:
FTic time; time.tic();
// check if aready set
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
rank = ComputeAndCompress(epsilon, U, C, B);
rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
......@@ -135,13 +135,13 @@ public:
* @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
* @param[out] B matrix of size \f$\ell^3\times r\f$
*/
static unsigned int ComputeAndCompress(const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
static unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
/**
* Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
* file
*/
static void ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon);
static void ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon);
/**
* Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
......@@ -182,19 +182,19 @@ public:
{
const unsigned int idx
= (transfer[0]+3)*7*7 + (transfer[1]+3)*7 + (transfer[2]+3);
const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(const unsigned int idx, FReal CellWidth,
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(FReal CellWidth,
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel.getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
}
......@@ -229,7 +229,8 @@ public:
template <int ORDER, class MatrixKernelClass>
unsigned int
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilon,
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const FReal epsilon,
FReal* &U,
FReal* &C,
FReal* &B)
......@@ -241,8 +242,6 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
FPoint X[nnodes], Y[nnodes];
// set roots of target cell (X)
FChebTensor<order>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
// init matrix kernel
const MatrixKernelClass MatrixKernel;
// allocate memory and compute 316 m2l operators
FReal *_U, *_C, *_B;
......@@ -260,7 +259,7 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
for (unsigned int n=0; n<nnodes; ++n)
for (unsigned int m=0; m<nnodes; ++m)
_C[counter*nnodes*nnodes + n*nnodes + m]
= MatrixKernel.evaluate(X[m], Y[n]);
= MatrixKernel->evaluate(X[m], Y[n]);
// increment interaction counter
counter++;
}
......@@ -330,14 +329,14 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
template <int ORDER, class MatrixKernelClass>
void
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const FReal epsilon)
FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon)
{
// measure time
FTic time; time.tic();
// start computing process
FReal *U, *C, *B;
U = C = B = nullptr;
const unsigned int rank = ComputeAndCompress(epsilon, U, C, B);
const unsigned int rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
// store into binary file
const std::string filename(getFileName(epsilon));
std::ofstream stream(filename.c_str(),
......
......@@ -56,6 +56,9 @@ class FChebSymKernel
typedef SymmetryHandler<ORDER, MatrixKernelClass::Type> SymmetryHandlerClass;
enum {nnodes = AbstractBaseClass::nnodes};
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for handling all symmetries
const FSmartPointer<SymmetryHandlerClass,FSmartPointerMemory> SymHandler;
......@@ -105,11 +108,13 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebSymKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal Epsilon)
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const MatrixKernelClass *const inMatrixKernel,
const FReal Epsilon)
: AbstractBaseClass(inTreeHeight, inBoxWidth, inBoxCenter),
SymHandler(new SymmetryHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon, inBoxWidth, inTreeHeight)),
MatrixKernel(inMatrixKernel),
SymHandler(new SymmetryHandlerClass(MatrixKernel, Epsilon, inBoxWidth, inTreeHeight)),
Loc(nullptr), Mul(nullptr), countExp(nullptr)
{
this->allocateMemoryForPermutedExpansions();
......@@ -126,8 +131,9 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebSymKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const MatrixKernelClass *const inMatrixKernel) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
{}
......@@ -135,6 +141,7 @@ public:
/** Copy constructor */
FChebSymKernel(const FChebSymKernel& other)
: AbstractBaseClass(other),
MatrixKernel(other.MatrixKernel),
SymHandler(other.SymHandler),
Loc(nullptr), Mul(nullptr), countExp(nullptr)
{
......@@ -264,7 +271,7 @@ public:
// multiply (mat-mat-mul)
FReal Compressed [nnodes * 24];
const FReal scale = AbstractBaseClass::MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
const FReal scale = MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
for (unsigned int pidx=0; pidx<343; ++pidx) {
const unsigned int count = countExp[pidx];
if (count) {
......@@ -462,14 +469,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
}
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,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
}
};
......
......@@ -303,7 +303,7 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal
FChebTensor<ORDER>::setRootOfWeights(weights);
// now the entry-computer is responsible for weighting the matrix entries
EntryComputer<MatrixKernelClass> Computer(nnodes, X, nnodes, Y, weights);
EntryComputer<MatrixKernelClass> Computer(MatrixKernel, nnodes, X, nnodes, Y, weights);
// start timer
time.tic();
......
......@@ -59,6 +59,9 @@ protected://PB: for OptiDis
typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for M2L operator
FSmartPointer< M2LHandlerClass,FSmartPointerMemory> M2LHandler;
......@@ -71,11 +74,12 @@ public:
FChebTensorialKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const MatrixKernelClass *const inMatrixKernel,
const FReal inBoxWidthExtension,
const FReal Epsilon,
const double inMatParam = 0.0)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension),
MatrixKernel(inMatrixKernel),
M2LHandler(new M2LHandlerClass(MatrixKernel,
inTreeHeight,
inBoxWidth,
inBoxWidthExtension,
......@@ -143,7 +147,7 @@ public:
// scale factor for homogeneous case
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(ExtendedCellWidth));
const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
......@@ -158,8 +162,7 @@ public:
const int idxMul = idxV*nRhs + idxRhs;
// get index in matrix kernel
const unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
const unsigned int d = MatrixKernel->getPosition(idxLhs);
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
......@@ -237,14 +240,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
}
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,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
}
};
......
......@@ -82,8 +82,6 @@ class FChebTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyabl
ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
ncmp = MatrixKernelClass::NCMP};
// const MatrixKernelClass MatrixKernel;
// FReal *U, *C, *B;
FReal *U, *B;
FReal** C;
......
......@@ -80,7 +80,11 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelR(const FReal = 0.0, const unsigned int = 0) {}
FInterpMatrixKernelR() {}
// copy ctor
FInterpMatrixKernelR(const FInterpMatrixKernelR& other)
{}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -142,9 +146,14 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
static const unsigned int NLHS = 1; //< dim of loc exp
FReal LX,LY,LZ ;
FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : LX(1.0),LY(1.0),LZ(1.0)
FInterpMatrixKernelRH() : LX(1.0),LY(1.0),LZ(1.0)
{ }
// copy ctor
FInterpMatrixKernelRH(const FInterpMatrixKernelRH& other)
: LX(other.LX), LY(other.LY), LZ(other.LZ)
{}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
......@@ -207,7 +216,11 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelRR(const FReal = 0.0, const unsigned int = 0) {}
FInterpMatrixKernelRR() {}
// copy ctor
FInterpMatrixKernelRR(const FInterpMatrixKernelRR& other)
{}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -272,7 +285,11 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FInterpMatrixKernelLJ(const FReal = 0.0, const unsigned int = 0) {}
FInterpMatrixKernelLJ() {}
// copy ctor
FInterpMatrixKernelLJ(const FInterpMatrixKernelLJ& other)
{}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -388,6 +405,11 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
: _i(indexTab[d]), _j(indexTab[d+NCMP]), _CoreWidth2(CoreWidth*CoreWidth)
{}
// copy ctor
FInterpMatrixKernel_R_IJ(const FInterpMatrixKernel_R_IJ& other)
: _i(other._i), _j(other._j), _CoreWidth2(other._CoreWidth2)
{}
// returns position in reduced storage from position in full 3x3 matrix
unsigned int getPosition(const unsigned int n) const
{return applyTab[n];}
......@@ -487,6 +509,11 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
: _i(indexTab[d]), _j(indexTab[d+NCMP]), _k(indexTab[d+2*NCMP]), _CoreWidth2(CoreWidth*CoreWidth)
{}
// copy ctor
FInterpMatrixKernel_R_IJK(const FInterpMatrixKernel_R_IJK& other)
: _i(other._i), _j(other._j), _k(other._k), _CoreWidth2(other._CoreWidth2)
{}
// returns position in reduced storage from position in full 3x3x3 matrix
unsigned int getPosition(const unsigned int n) const
{return applyTab[n];}
......@@ -598,10 +625,10 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
number of rows and cols and on the coordinates x and y and the type of the
generating matrix-kernel function.
*/
template <typename MatrixKernel>
template <typename MatrixKernelClass>
class EntryComputer
{
const MatrixKernel Kernel;
const MatrixKernelClass *const MatrixKernel;
const unsigned int nx, ny;
const FPoint *const px, *const py;
......@@ -609,12 +636,11 @@ class EntryComputer
const FReal *const weights;
public:
explicit EntryComputer(const unsigned int _nx, const FPoint *const _px,
explicit EntryComputer(const MatrixKernelClass *const inMatrixKernel,
const unsigned int _nx, const FPoint *const _px,
const unsigned int _ny, const FPoint *const _py,
const FReal *const _weights = NULL,
const unsigned int idxK = 0,
const FReal matparam = 0.0)
: Kernel(matparam,idxK), nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
const FReal *const _weights = NULL)
: MatrixKernel(inMatrixKernel), nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
void operator()(const unsigned int xbeg, const unsigned int xend,
const unsigned int ybeg, const unsigned int yend,
......@@ -624,11 +650,11 @@ public:
if (weights) {
for (unsigned int j=ybeg; j<yend; ++j)
for (unsigned int i=xbeg; i<xend; ++i)
data[idx++] = weights[i] * weights[j] * Kernel.evaluate(px[i], py[j]);
data[idx++] = weights[i] * weights[j] * MatrixKernel->evaluate(px[i], py[j]);
} else {
for (unsigned int j=ybeg; j<yend; ++j)
for (unsigned int i=xbeg; i<xend; ++i)
data[idx++] = Kernel.evaluate(px[i], py[j]);
data[idx++] = MatrixKernel->evaluate(px[i], py[j]);
}
/*
......
......@@ -49,8 +49,6 @@ protected:
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for P2P operator
const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Height of the entire oct-tree
const unsigned int TreeHeight;
/// Corner of oct-tree box
......@@ -62,9 +60,6 @@ protected:
/// Extension of the box width ( same for all level! )
const FReal BoxWidthExtension;
/// Parameter to pass to matrix kernel (material specific or anything)
const double MatParam;
/**
* Compute center of leaf cell from its tree coordinate.
* @param[in] Coordinate tree coordinate
......@@ -86,18 +81,15 @@ public:
FAbstractUnifKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal inBoxWidthExtension = 0.0,
const double inMatParam = 0.0)
const FReal inBoxWidthExtension = 0.0)
: Interpolator(new InterpolatorClass(inTreeHeight,
inBoxWidth,
inBoxWidthExtension)),
MatrixKernel(new MatrixKernelClass(inMatParam)),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
BoxWidthExtension(inBoxWidthExtension),
MatParam(inMatParam)
BoxWidthExtension(inBoxWidthExtension)
{
/* empty */
}
......
......@@ -51,6 +51,9 @@ class FUnifKernel
typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for M2L operator
const M2LHandlerClass M2LHandler;
......@@ -63,9 +66,11 @@ public:
*/
FUnifKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter)
const FPoint& inBoxCenter,
const MatrixKernelClass *const inMatrixKernel)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
MatrixKernel(inMatrixKernel),
M2LHandler(MatrixKernel,
inTreeHeight,
inBoxWidth)
{ }
......@@ -131,7 +136,7 @@ public:
const int TreeLevel)
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FComplexe *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);
......@@ -208,7 +213,7 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
}
......@@ -216,7 +221,7 @@ public:
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
}
};
......
......@@ -75,6 +75,9 @@ protected://PB: for OptiDis
typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for P2P and M2L operators
const MatrixKernelClass *const MatrixKernel;
/// Needed for M2L operator
const M2LHandlerClass M2LHandler;
......@@ -87,10 +90,11 @@ public:
FUnifTensorialKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal inBoxWidthExtension,
const double inMatParam = 0.0)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension,inMatParam),
M2LHandler(AbstractBaseClass::MatrixKernel.getPtr(),
const MatrixKernelClass *const inMatrixKernel,
const FReal inBoxWidthExtension)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inBoxWidthExtension),
MatrixKernel(inMatrixKernel),
M2LHandler(MatrixKernel,
inTreeHeight,
inBoxWidth,
inBoxWidthExtension)
......@@ -157,7 +161,7 @@ public:
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(ExtendedCellWidth));
const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
......@@ -175,8 +179,7 @@ public:
const int idxMul = idxV*nRhs + idxRhs;
// get index in matrix kernel
const unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxLhs);
const unsigned int d = MatrixKernel->getPosition(idxLhs);
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
......@@ -248,14 +251,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,MatrixKernel);
}
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,AbstractBaseClass::MatrixKernel.getPtr());
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,MatrixKernel);
}
};
......
......@@ -78,7 +78,7 @@ int main(int argc, char* argv[])
const unsigned int NLHS = MatrixKernelClass::NLHS;
const double CoreWidth = 0.1;
const MatrixKernelClass DirectMatrixKernel(CoreWidth);
const MatrixKernelClass MatrixKernel(CoreWidth);
std::cout<< "Interaction kernel: ";
if(MK_ID == ONE_OVER_R) std::cout<< "1/R";
else if(MK_ID == R_IJ)
......@@ -143,16 +143,16 @@ int main(int argc, char* argv[])
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
particles[idxOther].forces[0], particles[idxOther].forces[1],
particles[idxOther].forces[2], particles[idxOther].potential,
&DirectMatrixKernel);
else if(MK_ID == ONE_OVER_R)
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue[0],
particles[idxTarget].forces[0], particles[idxTarget].forces[1],
particles[idxTarget].forces[2], particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue[0],
particles[idxOther].forces[0], particles[idxOther].forces[1],
particles[idxOther].forces[2], particles[idxOther].potential);
&MatrixKernel);
// else if(MK_ID == ONE_OVER_R)
// FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
// particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue[0],
// particles[idxTarget].forces[0], particles[idxTarget].forces[1],
// particles[idxTarget].forces[2], particles[idxTarget].potential,
// particles[idxOther].position.getX(), particles[idxOther].position.getY(),
// particles[idxOther].position.getZ(), particles[idxOther].physicalValue[0],