Commit 1dde893e authored by berenger-bramas's avatar berenger-bramas
Browse files

Update the spherical class and the elec-foces-kernel

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@352 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 132cd57e
...@@ -111,7 +111,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta ...@@ -111,7 +111,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta
preM2MTransitionsPer = new FComplexe[periodicLevels * 8 * harmonic.getExpSize()]; preM2MTransitionsPer = new FComplexe[periodicLevels * 8 * harmonic.getExpSize()];
FReal treeWidthAtLevel = boxWidth; FReal treeWidthAtLevel = boxWidth;
for(int idxLevel = -1 ; idxLevel <= -periodicLevels ; --idxLevel ){ for(int idxLevel = -1 ; idxLevel >= -periodicLevels ; --idxLevel ){
const F3DPosition father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel); const F3DPosition father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel);
treeWidthAtLevel *= 2; treeWidthAtLevel *= 2;
...@@ -143,7 +143,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta ...@@ -143,7 +143,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta
// so 6 in each dimension // so 6 in each dimension
treeWidthAtLevel = boxWidth*2; treeWidthAtLevel = boxWidth*2;
preM2LTransitionsPer = new FComplexe[periodicLevels * (7 * 7 * 7) * devM2lP]; preM2LTransitionsPer = new FComplexe[periodicLevels * (7 * 7 * 7) * devM2lP];
for(int idxLevel = -1 ; idxLevel <= -periodicLevels ; --idxLevel ){ for(int idxLevel = -1 ; idxLevel >= -periodicLevels ; --idxLevel ){
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){ for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){ for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){ for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
...@@ -870,10 +870,11 @@ private: ...@@ -870,10 +870,11 @@ private:
// We want: - gradient(POTENTIAL_SIGN potential). // We want: - gradient(POTENTIAL_SIGN potential).
// The -(- 1.0) computing is not the most efficient programming ... // The -(- 1.0) computing is not the most efficient programming ...
const FReal signe = 1.0; const FReal signe = 1.0;
if( FMath::Epsilon < spherical.getR()){
force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR()); force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR());
force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR()); force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR());
force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta())); force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta()));
}
///////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////
//spherical_position_Set_ph //spherical_position_Set_ph
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define FSPHERICAL_HPP #define FSPHERICAL_HPP
#include "FGlobal.hpp" #include "FGlobal.hpp"
#include "FMath.hpp"
#include "F3DPosition.hpp" #include "F3DPosition.hpp"
/** This class is a Spherical position /** This class is a Spherical position
...@@ -20,9 +21,15 @@ public: ...@@ -20,9 +21,15 @@ public:
const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY()); const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY());
this->r = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ())); this->r = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ()));
this->phi = FMath::Atan2(inVector.getY(),inVector.getX()); this->phi = FMath::Atan2(inVector.getY(),inVector.getX());
if( r < FMath::Epsilon ){ // if r == 0 we cannot divide!
this->cosTheta = FReal(1);
this->sinTheta = FReal(1);
}
else {
this->cosTheta = inVector.getZ() / r; this->cosTheta = inVector.getZ() / r;
this->sinTheta = FMath::Sqrt(x2y2) / r; this->sinTheta = FMath::Sqrt(x2y2) / r;
} }
}
/** Get the rayon */ /** Get the rayon */
FReal getR() const{ FReal getR() const{
......
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