Commit b0252c07 authored by Matthias Messner's avatar Matthias Messner
parents 7eaa05a3 de73f612
This diff is collapsed.
......@@ -38,7 +38,7 @@ public:
* @param inImag the imaginary
* @param inReal the complex[0]
*/
explicit FComplexe(const FReal inImag, const FReal inReal) {
explicit FComplexe(const FReal inReal, const FReal inImag) {
complex[0] = inReal;
complex[1] = inImag;
}
......@@ -93,6 +93,16 @@ public:
this->complex[1] = inImag;
}
/** Return the conjugate */
FComplexe conjugate() const{
return FComplexe(complex[0],-complex[1]);
}
/** Return the conjugate */
FComplexe negate() const{
return FComplexe(-complex[0],-complex[1]);
}
/**
* Operator +=
* in complex[0] with other complex[0], same for complex[1]
......@@ -175,6 +185,12 @@ public:
this->complex[1] += (other.complex[0] * another.complex[1]) + (other.complex[1] * another.complex[0]);
}
/** Mul other and another and add the result to current complexe */
void equalMul(const FComplexe& other, const FComplexe& another){
this->complex[0] = (other.complex[0] * another.complex[0]) - (other.complex[1] * another.complex[1]);
this->complex[1] = (other.complex[0] * another.complex[1]) + (other.complex[1] * another.complex[0]);
}
/** To cast to FReal */
static FReal* ToFReal(FComplexe*const array){
return reinterpret_cast<FReal*>(array);
......
......@@ -80,8 +80,15 @@ struct FMath{
}
/** To get pow */
static double pow(double x, double y){
return ::pow(x,y);
}
static double pow(float x, float y){
return ::powf(x,y);
}
template <class NumType>
static NumType pow(const NumType inValue, int power){
if(power<0) power = -power;
NumType result = 1;
while(power-- > 0) result *= inValue;
return result;
......
......@@ -36,6 +36,11 @@ class FSpherical {
FReal theta; //!< the inclination angle
FReal phi; //!< the azimuth angle
public:
/** Default Constructor, set attributes to 0 */
FSpherical()
: r(0), cosTheta(0), sinTheta(0), theta(0), phi(0) {
}
/** From now, we just need a constructor based on a 3D position */
explicit FSpherical(const FPoint& inVector){
const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY());
......
// ===================================================================================
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
// ===================================================================================
#include <limits>
#include <iostream>
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Kernels/Spherical/FSphericalParticle.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FMemUtils.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
#include "../../Src/Core/FFmmAlgorithmTask.hpp"
#include "../../Src/Files/FFmaLoader.hpp"
/** We need to know the position of the particle in the array */
class IndexedParticle : public FSphericalParticle {
int index;
public:
IndexedParticle(): index(-1){}
int getIndex() const{
return index;
}
void setIndex( const int inIndex ){
index = inIndex;
}
};
int main(int argc, char** argv){
static const int P = 1;
typedef IndexedParticle ParticleClass;
typedef FRotationCell<P> CellClass;
typedef FVector<ParticleClass> ContainerClass;
typedef FSimpleLeaf<ParticleClass, ContainerClass > LeafClass;
typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FRotationKernel<ParticleClass, CellClass, ContainerClass , P> KernelClass;
typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassThread;
typedef FFmmAlgorithmTask<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassTask;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
std::cout << ">> You can pass -sequential or -task (thread by default).\n";
//////////////////////////////////////////////////////////////
const int NbLevels = FParameters::getValue(argc,argv,"-h", 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
FTic counter;
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
std::cout << "Opening : " << filename << "\n";
FFmaLoader<ParticleClass> loader(filename);
if(!loader.isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
// -----------------------------------------------------
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
// -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
//loader.fillTree(tree);
ParticleClass* const particles = new ParticleClass[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(particles[idxPart]);
particles[idxPart].setIndex( idxPart );
tree.insert(particles[idxPart]);
}
counter.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Create kernel ..." << std::endl;
counter.tic();
KernelClass kernels(P, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
counter.tac();
std::cout << "Done " << " in " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Working on particles ..." << std::endl;
if( FParameters::findParameter(argc,argv,"-sequential") != FParameters::NotFound){
FmmClass algo(&tree,&kernels);
counter.tic();
algo.execute();
}
else if( FParameters::findParameter(argc,argv,"-task") != FParameters::NotFound){
FmmClassTask algo(&tree,&kernels);
counter.tic();
algo.execute();
}
else {
FmmClassThread algo(&tree,&kernels);
counter.tic();
algo.execute();
}
counter.tac();
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
{ // get sum forces&potential
FReal potential = 0;
FPoint forces;
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
ContainerClass::ConstBasicIterator iter(*octreeIterator.getCurrentListTargets());
while( iter.hasNotFinished() ){
potential += iter.data().getPotential() * iter.data().getPhysicalValue();
forces += iter.data().getForces();
iter.gotoNext();
}
} while(octreeIterator.moveRight());
std::cout << "Foces Sum x = " << forces.getX() << " y = " << forces.getY() << " z = " << forces.getZ() << std::endl;
std::cout << "Potential = " << potential << std::endl;
}
// -----------------------------------------------------
{
std::cout << "Compute direct interaction for all\n";
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
kernels.particlesMutualInteraction(&particles[idxTarget], &particles[idxOther]);
}
}
std::cout << "Compute Diff...\n";
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
typename ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());
while( leafIter.hasNotFinished() ){
const ParticleClass& other = particles[leafIter.data().getIndex()];
potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
std::cout << "Direct Potential = " << other.getPotential() << " Fmm Potential = " << leafIter.data().getPotential() << std::endl; // Remove Me
fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
fz.add(other.getForces().getZ(),leafIter.data().getForces().getZ());
leafIter.gotoNext();
}
} while(octreeIterator.moveRight());
}
// Print for information
std::cout << "Potential diff is = " << potentialDiff.getL2Norm() << " " << potentialDiff.getInfNorm() << std::endl;
std::cout << "Fx diff is = " << fx.getL2Norm() << " " << fx.getInfNorm() << std::endl;
std::cout << "Fy diff is = " << fy.getL2Norm() << " " << fy.getInfNorm() << std::endl;
std::cout << "Fz diff is = " << fz.getL2Norm() << " " << fz.getInfNorm() << std::endl;
}
delete[] particles;
return 0;
}
This diff is collapsed.
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment