Commit f81a1e74 authored by COULAUD Olivier's avatar COULAUD Olivier

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

# By Pierre Blanchard
# Via Pierre Blanchard
* 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm:
  Provided generic functions for the evaluation of the direct interactions (pass matrix kernel as argument).
parents 0deb8d64 728d90e3
......@@ -29,6 +29,7 @@
#include "Files/FFmaGenericLoader.hpp"
#include "Kernels/P2P/FP2P.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
//
......@@ -118,6 +119,9 @@ int main(int argc, char ** argv){
// ----------------------------------------------------------------------------------------------------------
// COMPUTATION
// ----------------------------------------------------------------------------------------------------------
// interaction kernel evaluator
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
FReal denergy = 0.0;
//
// computation
......@@ -142,7 +146,7 @@ int main(int argc, char ** argv){
particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
particles[idxTarget].getPosition().getZ(),particles[idxTarget].getPhysicalValue(),
&particles[idxTarget].setForces()[0],&particles[idxTarget].setForces()[1],
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential());
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),&MatrixKernel);
}
}
} // end for
......
......@@ -212,14 +212,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());
}
};
......
......@@ -462,14 +462,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,6 +86,7 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -94,11 +95,29 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
block[0] = one_over_r;
blockDerivative[0][0] = one_over_r3 * xy.getX();
blockDerivative[0][1] = one_over_r3 * xy.getY();
blockDerivative[0][2] = one_over_r3 * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -126,6 +145,7 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : LX(1.0),LY(1.0),LZ(1.0)
{ }
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -145,6 +165,23 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal one_over_rL = FReal(1.) / FMath::Sqrt(LX*xy.getX()*xy.getX() +
LY*xy.getY()*xy.getY() +
LZ*xy.getZ()*xy.getZ());
const FReal one_over_rL3 = one_over_rL*one_over_rL*one_over_rL;
block[0] = one_over_rL;
blockDerivative[0][0] = FMath::Sqrt(LX) * one_over_rL3 * xy.getX();
blockDerivative[0][1] = FMath::Sqrt(LY) * one_over_rL3 * xy.getY();
blockDerivative[0][2] = FMath::Sqrt(LZ) * one_over_rL3 * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -176,6 +213,7 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -184,11 +222,31 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal r2 = FReal(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
const FReal one_over_r2 = FReal(1.) / (r2);
const FReal one_over_r4 = one_over_r2*one_over_r2;
block[0] = one_over_r2;
const FReal coef = FReal(-2.) * one_over_r4;
blockDerivative[0][0] = coef * xy.getX();
blockDerivative[0][1] = coef * xy.getY();
blockDerivative[0][2] = coef * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -220,6 +278,7 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -231,11 +290,31 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
return one_over_r6 * one_over_r6 - one_over_r6;
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal r = xy.norm();
const FReal r2 = r*r;
const FReal r3 = r2*r;
const FReal one_over_r6 = FReal(1.) / (r3*r3);
const FReal one_over_r8 = one_over_r6 / (r2);
block[0] = one_over_r6 * one_over_r6 - one_over_r6;
const FReal coef = FReal(12.0)*one_over_r6*one_over_r8 - FReal(6.0)*one_over_r8;
blockDerivative[0][0]= coef * xy.getX();
blockDerivative[0][1]= coef * xy.getY();
blockDerivative[0][2]= coef * xy.getZ();
}
FReal getScaleFactor(const FReal, const int) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
......@@ -317,6 +396,7 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
FReal getCoreWidth2() const
{return _CoreWidth2;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -343,6 +423,7 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
const FPoint xy(x-y);
......@@ -414,6 +495,7 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
FReal getCoreWidth2() const
{return _CoreWidth2;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
// Convention for anti-symmetric kernels xy=y-x instead of xy=x-y
......@@ -460,6 +542,7 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
// Convention for anti-symmetric kernels xy=y-x instead of xy=x-y
......
......@@ -4,69 +4,51 @@
#include "../P2P/FP2P.hpp"
template <KERNEL_FUNCTION_IDENTIFIER Identifier, int NVALS>
struct DirectInteractionComputer;
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteractionComputer<ONE_OVER_R, 1>
template <KERNEL_FUNCTION_IDENTIFIER Identifier, int NVALS>
struct DirectInteractionComputer
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const MatrixKernel){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14,MatrixKernel);
}
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
const int inSize,
const MatrixKernelClass *const MatrixKernel){
FP2P::FullRemote(inTargets,inNeighbors,inSize,MatrixKernel);
}
};
//! Specialization for Laplace potential on Non uniform domain
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteractionComputer<ONE_OVER_RH, 1>
struct DirectInteractionComputer<ONE_OVER_R, 1>
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const /*MatrixKernel*/){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
//
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
const int inSize,
const MatrixKernelClass *const /*MatrixKernel*/){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for Lennard-Jones potential */
template <>
struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, 1>
{
template <typename ContainerClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote(ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for GradGradR potential */
template <>
struct DirectInteractionComputer<R_IJ, 1>
......@@ -116,18 +98,20 @@ struct DirectInteractionComputer<R_IJK, 1>
template <int NVALS>
struct DirectInteractionComputer<ONE_OVER_R, NVALS>
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const /*MatrixKernel*/){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote(ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
const int inSize,
const MatrixKernelClass *const /*MatrixKernel*/){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
......@@ -135,29 +119,6 @@ struct DirectInteractionComputer<ONE_OVER_R, NVALS>
};
/*! Specialization for Lennard-Jones potential */
template <int NVALS>
struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, NVALS>
{
template <typename ContainerClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualLJ(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::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for GradGradR potential */
template <int NVALS>
struct DirectInteractionComputer<R_IJ, NVALS>
......
This diff is collapsed.
......@@ -208,14 +208,15 @@ 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);
ContainerClass* const inNeighbors[27], const int /*inSize*/)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
}
};
......
......@@ -78,6 +78,7 @@ int main(int argc, char ** argv){
const unsigned int ORDER = 12;
// typedefs
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
......@@ -236,7 +237,7 @@ int main(int argc, char ** argv){
particlesDirect[idxTarget].position.getX(), particlesDirect[idxTarget].position.getY(),
particlesDirect[idxTarget].position.getZ(),particlesDirect[idxTarget].physicalValue,
&particlesDirect[idxTarget].forces[0],&particlesDirect[idxTarget].forces[1],
&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential);
&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential, &MatrixKernel);
}
}
//
......@@ -258,7 +259,7 @@ int main(int argc, char ** argv){
FP2P::NonMutualParticles(
source.position.getX(), source.position.getY(),source.position.getZ(),source.physicalValue,
particlesDirect[idxTarget].position.getX(), particlesDirect[idxTarget].position.getY(),particlesDirect[idxTarget].position.getZ(),particlesDirect[idxTarget].physicalValue,
&particlesDirect[idxTarget].forces[0],&particlesDirect[idxTarget].forces[1],&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential
&particlesDirect[idxTarget].forces[0],&particlesDirect[idxTarget].forces[1],&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential, &MatrixKernel
);
}
}
......
......@@ -91,6 +91,7 @@ void doATest(const int NbParticles, const int minP, const int maxP, const int mi
if(computeDirectAndDiff){
printf("Compute direct!\n");
counter.tic();
const FInterpMatrixKernelR MatrixKernel;
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
......@@ -100,7 +101,7 @@ void doATest(const int NbParticles, const int minP, const int maxP, const int mi
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
&particles[idxOther].forces[0],&particles[idxOther].forces[1],
&particles[idxOther].forces[2],&particles[idxOther].potential);
&particles[idxOther].forces[2],&particles[idxOther].potential,&MatrixKernel);
}
}
if(timeForDirect) *timeForDirect = counter.tacAndElapsed();
......@@ -299,6 +300,7 @@ int main(int argc, char ** argv){
FmmClass algo(&tree,&kernels);
algo.execute();
const FInterpMatrixKernelR MatrixKernel;
FP2P::MutualParticles(centeredParticle.position.getX(), centeredParticle.position.getY(),
centeredParticle.position.getZ(),centeredParticle.physicalValue,
&centeredParticle.forces[0],&centeredParticle.forces[1],
......@@ -306,7 +308,7 @@ int main(int argc, char ** argv){
otherParticle.position.getX(), otherParticle.position.getY(),
otherParticle.position.getZ(),otherParticle.physicalValue,
&otherParticle.forces[0],&otherParticle.forces[1],
&otherParticle.forces[2],&otherParticle.potential);
&otherParticle.forces[2],&otherParticle.potential,&MatrixKernel);
{ // Check that each particle has been summed with all other
......
......@@ -279,6 +279,8 @@ int main(int argc, char ** argv){
particlesDirect = new EwalParticle[loader->getNumberOfParticles()];
const FInterpMatrixKernelR MatrixKernel;
FReal denergy = 0.0;
FMath::FAccurater dfx, dfy, dfz ;
//
......@@ -305,7 +307,7 @@ int main(int argc, char ** argv){
part.position.getX(), part.position.getY(),
part.position.getZ(),part.physicalValue,
&part.forces[0],&part.forces[1],
&part.forces[2],&part.potential);
&part.forces[2],&part.potential,&MatrixKernel);
}
}
//
......@@ -333,7 +335,7 @@ int main(int argc, char ** argv){
FP2P::NonMutualParticles(
source.position.getX(), source.position.getY(),source.position.getZ(),source.physicalValue,
part.position.getX(), part.position.getY(),part.position.getZ(),part.physicalValue,
&part.forces[0],&part.forces[1],&part.forces[2],&part.potential
&part.forces[0],&part.forces[1],&part.forces[2],&part.potential,&MatrixKernel
);
}
// std::cout <<std::endl<<std::endl<<std::endl;
......
......@@ -90,6 +90,9 @@ int main(int argc, char* argv[])
particles[idxPart].forces[2] = 0.0;
}
{
// interaction evaluator
const FInterpMatrixKernelR MatrixKernel;
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
......@@ -99,7 +102,7 @@ int main(int argc, char* argv[])
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
&particles[idxOther].forces[0],&particles[idxOther].forces[1],
&particles[idxOther].forces[2],&particles[idxOther].potential);
&particles[idxOther].forces[2],&particles[idxOther].potential, &MatrixKernel);
}
}
}
......
......@@ -69,6 +69,12 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// interaction kernel evaluator
// typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FInterpMatrixKernelRR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
// init particles position and physical value
struct TestParticle{
FPoint position;
......@@ -110,7 +116,8 @@ int main(int argc, char* argv[])
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential);
&particles[idxOther].forces[2], &particles[idxOther].potential,
&MatrixKernel);
}
}
}
......@@ -130,8 +137,6 @@ int main(int argc, char* argv[])
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
// typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FUnifCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
......
......@@ -211,7 +211,10 @@ int main(int argc, char* argv[])
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
if(NPV==3) // R_IJ or IOR
if(NPV==1) // scalar kernels like ONE_OVER_R
tree.insert(particles[idxPart].position, idxPart,
particles[idxPart].physicalValue[0]);
else if(NPV==3) // R_IJ or IOR
tree.insert(particles[idxPart].position, idxPart,
particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
else if(NPV==9) // R_IJK
......
......@@ -64,6 +64,11 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// interaction kernel evaluator
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
// init particles position and physical value
struct TestParticle{
FPoint position;
......@@ -105,7 +110,7 @@ int main(int argc, char* argv[])
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential);
&particles[idxOther].forces[2], &particles[idxOther].potential,&MatrixKernel);
}
}
}
......@@ -124,8 +129,6 @@ int main(int argc, char* argv[])
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
......
......@@ -73,6 +73,10 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// interaction kernel evaluator
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
struct TestParticle{
FPoint position;
FReal forces[3];
......@@ -106,7 +110,7 @@ int main(int argc, char* argv[])
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
&particles[idxOther].forces[0],&particles[idxOther].forces[1],
&particles[idxOther].forces[2],&particles[idxOther].potential);
&particles[idxOther].forces[2],&particles[idxOther].potential,&MatrixKernel);
}
}
}
......@@ -121,7 +125,6 @@ int main(int argc, char* argv[])
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
......
......@@ -127,6 +127,7 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
// Run direct computation
//Print("Direct...");
const FInterpMatrixKernelR MatrixKernel;
for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < nbParticles ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
......@@ -136,7 +137,7 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
&particles[idxOther].forces[0],&particles[idxOther].forces[1],
&particles[idxOther].forces[2],&particles[idxOther].potential);
&particles[idxOther].forces[2],&particles[idxOther].potential,&MatrixKernel);
}
}
......
......@@ -18,26 +18,26 @@
// @FUSE_BLAS
// ================
#include "../Src/Utils/FGlobal.hpp"
#include "Utils/FGlobal.hpp"
#include "../Src/Containers/FOctree.hpp"
#include "../Src/Containers/FVector.hpp"
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"