Commit 596ca731 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
parents 4c09a276 ca0bf8b8
This diff is collapsed.
......@@ -18,48 +18,49 @@
#include "../../Components/FBasicParticleContainer.hpp"
class FP2PParticleContainer : public FBasicParticleContainer<5> {
typedef FBasicParticleContainer<5> Parent;
template<int NRHS = 1, int NLHS = 1>
class FP2PParticleContainer : public FBasicParticleContainer<NRHS+4*NLHS> {
typedef FBasicParticleContainer<NRHS+4*NLHS> Parent;
public:
FReal* getPhysicalValues(){
return Parent::getAttribute<0>();
FReal* getPhysicalValues(const int idxRhs = 0){
return Parent::getAttribute(0+idxRhs);
}
const FReal* getPhysicalValues() const {
return Parent::getAttribute<0>();
const FReal* getPhysicalValues(const int idxRhs = 0) const {
return Parent::getAttribute(0+idxRhs);
}
FReal* getPotentials(){
return Parent::getAttribute<1>();
FReal* getPotentials(const int idxLhs = 0){
return Parent::getAttribute(NRHS+idxLhs);
}
const FReal* getPotentials() const {
return Parent::getAttribute<1>();
const FReal* getPotentials(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+idxLhs);
}
FReal* getForcesX(){
return Parent::getAttribute<2>();
FReal* getForcesX(const int idxLhs = 0){
return Parent::getAttribute(NRHS+NLHS+idxLhs);
}
const FReal* getForcesX() const {
return Parent::getAttribute<2>();
const FReal* getForcesX(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+NLHS+idxLhs);
}
FReal* getForcesY(){
return Parent::getAttribute<3>();
FReal* getForcesY(const int idxLhs = 0){
return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
}
const FReal* getForcesY() const {
return Parent::getAttribute<3>();
const FReal* getForcesY(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
}
FReal* getForcesZ(){
return Parent::getAttribute<4>();
FReal* getForcesZ(const int idxLhs = 0){
return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
}
const FReal* getForcesZ() const {
return Parent::getAttribute<4>();
const FReal* getForcesZ(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
}
};
......
......@@ -6,8 +6,9 @@
#include "FP2PParticleContainer.hpp"
#include "../../Components/FParticleType.hpp"
class FP2PParticleContainerIndexed : public FP2PParticleContainer {
typedef FP2PParticleContainer Parent;
template<int NRHS = 1, int NLHS = 1>
class FP2PParticleContainerIndexed : public FP2PParticleContainer<NRHS,NLHS> {
typedef FP2PParticleContainer<NRHS,NLHS> Parent;
FVector<int> indexes;
......
......@@ -362,14 +362,13 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
const ContainerClass *const inParticles) const
{
// set all multipole expansions to zero
FBlas::setzero(/*nRhs**/nnodes, multipoleExpansion);
FBlas::setzero(nRhs*nnodes, multipoleExpansion);
// allocate stuff
const map_glob_loc map(center, width);
FPoint localPosition;
// loop over source particles
const FReal*const physicalValues = inParticles->getPhysicalValues(); // PB: TODO add multidim PhysVal
const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2];
......@@ -385,24 +384,24 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops
}
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M
// // this avoids nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// read physicalValue
const FReal weight = physicalValues[idxPart/*+idxRhs*nParticles*/]; //PB: TODO select next compo
const FReal*const physicalValues = inParticles->getPhysicalValues(idxRhs);
// compute weight
const FReal weight = physicalValues[idxPart];
// assemble multipole expansions
for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int idx = /*idxRhs*nnodes +*/ k*ORDER*ORDER + j*ORDER + i;
const unsigned int idx = idxRhs*nnodes + k*ORDER*ORDER + j*ORDER + i;
multipoleExpansion[idx] += L_of_x[i][0] * L_of_x[j][1] * L_of_x[k][2] * weight; // 3 * ORDER*ORDER*ORDER flops
}
}
}
// } // idxRhs
} // idxRhs
} // flops: N * (3 * ORDER*ORDER*ORDER + 3 * 3 * ORDER*(ORDER-1)) flops
......@@ -423,11 +422,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
const map_glob_loc map(center, width);
FPoint localPosition;
//const FReal*const physicalValues = inParticles->getPhysicalValues();
const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2];
FReal*const potentials = inParticles->getPotentials();
const unsigned int nParticles = inParticles->getNbParticles();
......@@ -444,9 +441,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops
}
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M
// // this avoid nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// interpolate and increment target value
FReal targetValue=0.;
......@@ -454,7 +449,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
for (unsigned int l=0; l<ORDER; ++l) {
for (unsigned int m=0; m<ORDER; ++m) {
for (unsigned int n=0; n<ORDER; ++n) {
const unsigned int idx = /*idxLhs*nnodes +*/ n*ORDER*ORDER + m*ORDER + l;
const unsigned int idx = idxLhs*nnodes + n*ORDER*ORDER + m*ORDER + l;
targetValue +=
L_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx];
} // ORDER * 4 flops
......@@ -462,10 +457,12 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
} // ORDER * ORDER * ORDER * 4 flops
}
// set potential
potentials[idxPart/*+idxLhs*nParticles*/] += (targetValue);
// get potential
FReal*const potentials = inParticles->getPotentials(idxLhs);
// add contribution to potential
potentials[idxPart] += (targetValue);
// } // idxLhs
} // idxLhs
} // N * (4 * ORDER * ORDER * ORDER + 9 * ORDER*(ORDER-1) ) flops
}
......@@ -495,14 +492,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
FReal L_of_x[ORDER][3];
FReal dL_of_x[ORDER][3];
const FReal*const physicalValues = inParticles->getPhysicalValues();
const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2];
FReal*const forcesX = inParticles->getForcesX();
FReal*const forcesY = inParticles->getForcesY();
FReal*const forcesZ = inParticles->getForcesZ();
//FReal*const potentials = inParticles->getPotentials();
// const unsigned int nParticles = inParticles->getNbParticles();
......@@ -521,9 +513,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
dL_of_x[o][2] = BasisType::dL(o, localPosition.getZ()); // TODO verify 3 * ORDER*(ORDER-1) flops
}
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M
// // this avoid nLhs/nLhs evaluations of the same interpolating polynomials
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// interpolate and increment forces value
FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)};
......@@ -531,7 +521,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
for (unsigned int l=0; l<ORDER; ++l) {
for (unsigned int m=0; m<ORDER; ++m) {
for (unsigned int n=0; n<ORDER; ++n) {
const unsigned int idx = /*idxLhs*nnodes +*/ n*ORDER*ORDER + m*ORDER + l;
const unsigned int idx = idxLhs*nnodes + n*ORDER*ORDER + m*ORDER + l;
forces[0] +=
dL_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx];
forces[1] +=
......@@ -545,19 +535,23 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
// scale forces
forces[0] *= jacobian[0];// / nnodes;
forces[1] *= jacobian[1];// / nnodes;
forces[2] *= jacobian[2];// / nnodes;
forces[0] *= jacobian[0];
forces[1] *= jacobian[1];
forces[2] *= jacobian[2];
}
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);
// set computed forces
// const unsigned int idx = idxPart+idxLhs*nParticles
forcesX[idxPart] += forces[0] * physicalValues[idxPart]; //PB: TODO update next compo
forcesY[idxPart] += forces[1] * physicalValues[idxPart]; //PB: TODO update next compo
forcesZ[idxPart] += forces[2] * physicalValues[idxPart]; //PB: TODO update next compo
}
forcesX[idxPart] += forces[0] * physicalValues[idxPart];
forcesY[idxPart] += forces[1] * physicalValues[idxPart];
forcesZ[idxPart] += forces[2] * physicalValues[idxPart];
}// idxLhs
// }// idxLhs
}
}
......
......@@ -41,11 +41,8 @@ class FTreeCoordinate;
* 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs (NVALS)
* are considered 2 separate features and are currently combined.
*
* 2) The present tensorial version is the most naive one. All tensorial aspects
* are handled in the kernel. A more optimal version would be to consider looping
* over nRhs/nLhs inside Interpolator::applyP2M/L2P in order to avoid extra
* evaluation of the interpolating polynomials. When it comes to applying M2L it is
* NOT much faster to loop over NRHSxNLHS inside applyM2L (at least for the Lagrange case)
* 2) When it comes to applying M2L it is NOT much faster to loop over NRHSxNLHS
* inside applyM2L (at least for the Lagrange case).
* 2-bis) The evaluation of the kernel matrix (see M2LHandler) should be done at once
* instead of compo-by-compo (TODO). On the other hand, the ChebyshevSym tensorial kernel
* requires the matrix kernel to be evaluated compo-by-compo since we currently use a scalar ACA.
......@@ -102,15 +99,17 @@ public:
const ContainerClass* const SourceParticles)
{
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxV*nRhs), SourceParticles);
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// update multipole index
int idxMul = idxV*nRhs + idxRhs;
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxMul), SourceParticles);
// 2) apply Discrete Fourier Transform
M2LHandler->applyZeroPaddingAndDFT(LeafCell->getMultipole(idxMul),
LeafCell->getTransformedMultipole(idxMul));
......@@ -211,15 +210,16 @@ public:
M2LHandler->unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxLoc),
const_cast<CellClass*>(LeafCell)->getLocal(idxLoc));
// 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxLoc), TargetParticles);
}
// 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxV*nLhs), TargetParticles);
// 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxLoc), TargetParticles);
// 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxV*nLhs), TargetParticles);
}
}// NVALS
}
......
......@@ -295,7 +295,7 @@ public:
int main(int argc, char* argv[]){
static const int P = 9;
typedef FRotationCell<P> CellClass;
typedef FP2PParticleContainer ContainerClass;
typedef FP2PParticleContainer<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
......
......@@ -124,7 +124,7 @@ int main(int argc, char* argv[])
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
// typedefs
typedef FP2PParticleContainerIndexed ContainerClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
......
......@@ -58,7 +58,7 @@ int main(int argc, char* argv[])
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
typedef FP2PParticleContainerIndexed ContainerClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
......
......@@ -66,12 +66,23 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// typedefs
// typedef FInterpMatrixKernelLJ MatrixKernelClass;
// typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FInterpMatrixKernel_R_IJ MatrixKernelClass; // not working with Non-Symmetric variant ! Because of UCB decomposition.
// and not working either with Symmetric variant because not symmetric...
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
// const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
// init particles position and physical value
struct TestParticle{
FPoint position;
FReal forces[3];
FReal physicalValue;
FReal potential;
FReal forces[3][NLHS];
FReal physicalValue[NRHS];
FReal potential[NLHS];
};
// open particle file
......@@ -85,11 +96,14 @@ int main(int argc, char* argv[])
loader.fillParticle(&position,&physicalValue);
// get copy
particles[idxPart].position = position;
particles[idxPart].physicalValue = physicalValue;
particles[idxPart].potential = 0.0;
particles[idxPart].forces[0] = 0.0;
particles[idxPart].forces[1] = 0.0;
particles[idxPart].forces[2] = 0.0;
for(unsigned idxRhs = 0; idxRhs<NRHS;++idxRhs)
particles[idxPart].physicalValue[idxRhs] = physicalValue;
for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
particles[idxPart].potential[idxLhs] = 0.0;
particles[idxPart].forces[0][idxLhs] = 0.0;
particles[idxPart].forces[1][idxLhs] = 0.0;
particles[idxPart].forces[2][idxLhs] = 0.0;
}
}
////////////////////////////////////////////////////////////////////
......@@ -110,12 +124,12 @@ int main(int argc, char* argv[])
// &particles[idxOther].forces[2], &particles[idxOther].potential);
FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
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,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential);
particles[idxOther].forces[0], particles[idxOther].forces[1],
particles[idxOther].forces[2], particles[idxOther].potential);
}
}
}
......@@ -134,18 +148,7 @@ int main(int argc, char* argv[])
const FReal epsilon = FReal(1e-8);
// typedefs
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
// typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FInterpMatrixKernel_R_IJ MatrixKernelClass; // not working with Non-Symmetric variant ! Because of UCB decomposition.
// and not working either with Symmetric variant because of lack of symmetry...
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
typedef FP2PParticleContainerIndexed ContainerClass;
// const unsigned int NDIM = NRHS + 4*NLHS;
// typedef FP2PTensorialParticleContainerIndexed<NDIM> ContainerClass; // TODO fix a TensorialParticleContainer for easy access to multidim PhysVal
typedef FP2PParticleContainerIndexed<NRHS,NLHS> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FChebCell<ORDER,NRHS,NLHS> CellClass;
......@@ -166,7 +169,8 @@ int main(int argc, char* argv[])
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
// PB: here we have to know the number of NRHS...
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
}
time.tac();
......@@ -187,44 +191,68 @@ int main(int argc, char* argv[])
{ // -----------------------------------------------------
std::cout << "\nError computation ... " << std::endl;
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
FMath::FAccurater potentialDiff[NLHS];
FMath::FAccurater fx[NLHS], fy[NLHS], fz[NLHS];
FReal checkPotential[20000];
FReal checkPhysVal[20000][NRHS];
FReal checkPotential[20000][NLHS];
FReal checkfx[20000][NLHS];
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int indexPartOrig = indexes[idxPart];
//PB: store potential in nbParticles array
checkPotential[indexPartOrig]=potentials[idxPart];
potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
const FReal*const physVals = leaf->getTargets()->getPhysicalValues(idxLhs);
const FReal*const potentials = leaf->getTargets()->getPotentials(idxLhs);
const FReal*const forcesX = leaf->getTargets()->getForcesX(idxLhs);
const FReal*const forcesY = leaf->getTargets()->getForcesY(idxLhs);
const FReal*const forcesZ = leaf->getTargets()->getForcesZ(idxLhs);
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int indexPartOrig = indexes[idxPart];
//PB: store potential in array[nbParticles]
checkPhysVal[indexPartOrig][idxLhs]=physVals[idxPart];
checkPotential[indexPartOrig][idxLhs]=potentials[idxPart];
checkfx[indexPartOrig][idxLhs]=forcesX[idxPart];
potentialDiff[idxLhs].add(particles[indexPartOrig].potential[idxLhs],potentials[idxPart]);
fx[idxLhs].add(particles[indexPartOrig].forces[0][idxLhs],forcesX[idxPart]);
fy[idxLhs].add(particles[indexPartOrig].forces[1][idxLhs],forcesY[idxPart]);
fz[idxLhs].add(particles[indexPartOrig].forces[2][idxLhs],forcesZ[idxPart]);
}
}
});
}
// std::cout << "Check Potential " << std::endl;
// std::cout << "Check Potential, forceX " << std::endl;
// for(int idxPart = 0 ; idxPart < 20 ; ++idxPart)
// std::cout << checkPotential[idxPart] << ", "<< particles[idxPart].potential << std::endl;
// for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
// std::cout << checkPhysVal[idxPart][idxLhs] << ", "<< particles[idxPart].physicalValue[idxLhs]<< "|| ";
// std::cout << checkPotential[idxPart][idxLhs] << ", "<< particles[idxPart].potential[idxLhs]<< "|| ";
// std::cout << checkfx[idxPart][idxLhs] << ", "<< particles[idxPart].forces[0][idxLhs] << std::endl;
// }
// std::cout << std::endl;
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
std::cout << "\nTest Chebyshev Tensorial returns wrong results! ChebInterpolator still needs to be modified like UnifInterpolator." << std::endl;
// std::cout << "\nAbsolute errors: " << std::endl;
// std::cout << "Potential: ";
// for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << potentialDiff[idxLhs] << ", " ;
// std::cout << std::endl;
// std::cout << "Fx: ";
// for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fx[idxLhs] << ", " ;
// std::cout << std::endl;
// std::cout << "Fy: ";
// for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fy[idxLhs] << ", " ;
// std::cout << std::endl;
// std::cout << "Fz: ";
// for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs) std::cout << fz[idxLhs] << ", " ;
// std::cout << std::endl;
} // -----------------------------------------------------
} // end Lagrange kernel
......
......@@ -119,7 +119,7 @@ int main(int argc, char* argv[])
const FReal epsilon = FReal(1e-7);
// typedefs
typedef FP2PParticleContainerIndexed ContainerClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
......@@ -199,7 +199,7 @@ int main(int argc, char* argv[])
// typedefs
typedef FSphericalCell CellClass;
typedef FP2PParticleContainerIndexed ContainerClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
//typedef FSphericalBlasKernel< CellClass, ContainerClass > KernelClass;
......
......@@ -55,9 +55,9 @@ int main(int argc, char* argv[])
// init timer
FTic time;
// typedefs
typedef FP2PParticleContainer ContainerClass;
typedef FP2PParticleContainer<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
......
......@@ -138,21 +138,21 @@ int main(int argc, char* argv[])
////////////////////////////////////////////////////////////////////
//
#ifdef ScalFMM_USE_BLAS
{ // begin Chebyshef kernel
{ // begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
// typedefs
typedef FP2PParticleContainerIndexed ContainerClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
// typedef FChebKernDel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
// typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
......@@ -221,7 +221,7 @@ int main(int argc, char* argv[])
// typedefs