Commit f9ebf2ba authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

kernel corrected

parent 942183b8
...@@ -36,9 +36,9 @@ protected: ...@@ -36,9 +36,9 @@ protected:
static const int LocalSize = ((P+1)*(P+2)*(P+3))/6; static const int LocalSize = ((P+1)*(P+2)*(P+3))/6;
//Multipole vector //Multipole vector
FVector<FReal> multipole_exp = FVector(MultipoleSize); FVector<FReal> multipole_exp = FVector<FReal>(MultipoleSize);
//Local vector //Local vector
FVector<FReal> local_exp = FVector(LocalSize); FVector<FReal> local_exp = FVector<FReal>(LocalSize);
public: public:
/** /**
...@@ -48,13 +48,13 @@ public: ...@@ -48,13 +48,13 @@ public:
} }
//Get multipole Vector //Get multipole Vector
FVector * getMultipole(void) FVector<FReal> * getMultipole(void)
{ {
return multipole_exp; return multipole_exp;
} }
//Get local Vector //Get local Vector
FVector * getLocal(void) FVector<FReal> * getLocal(void)
{ {
return local_exp; return local_exp;
} }
......
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
#define FTAYLORKERNEL_HPP #define FTAYLORKERNEL_HPP
#include "../../Components/FAbstractKernels.hpp" #include "../../Components/FAbstractKernels.hpp"
#include "../../Utils/FMemUtils.hpp"
/** /**
* @author Cyrille Piacibello * @author Cyrille Piacibello
...@@ -29,8 +29,8 @@ ...@@ -29,8 +29,8 @@
*/ */
template<class ParticleClass, class CellClass, class ContainerClass, int P> template< class ParticleClass, class CellClass, class ContainerClass, int P>
class FTaylorKernel : public FAbstractKernel<ParticleClass,CellClass,ContainerClass>{ class FTaylorKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> {
private: private:
//Size of the multipole and local vectors //Size of the multipole and local vectors
...@@ -62,7 +62,7 @@ private: ...@@ -62,7 +62,7 @@ private:
* Result : ...,[2,0,0],[1,1,0],[1,0,1]... 3-tuple are sorted * Result : ...,[2,0,0],[1,1,0],[1,0,1]... 3-tuple are sorted
* by size and alphabetical order. * by size and alphabetical order.
*/ */
void incPowers(int * const restrict a, int *const restrict b, int *const restict c) void incPowers(int * const FRestrict a, int *const FRestrict b, int *const FRestrict c)
{ {
int t = (*a)+(*b)+(*c); int t = (*a)+(*b)+(*c);
if(t==0) if(t==0)
...@@ -144,9 +144,9 @@ private: ...@@ -144,9 +144,9 @@ private:
tab[0]=1/(FMath::Sqrt(dx*dx+dy*dy+dz*dz)); tab[0]=1/(FMath::Sqrt(dx*dx+dy*dy+dz*dz));
FReal R3 = tab[0]/(dx*dx+dy*dy+dz*dz); FReal R3 = tab[0]/(dx*dx+dy*dy+dz*dz);
tab[1]=temp*dx; tab[1]=dx*R3;
tab[2]=temp*dy; tab[2]=dy*R3;
tab[3]=temp*dz; tab[3]=dz*R3;
FReal R5 = R3/(dx*dx+dy*dy+dz*dz); FReal R5 = R3/(dx*dx+dy*dy+dz*dz);
tab[4] = R3-dx*dx*R5; tab[4] = R3-dx*dx*R5;
tab[5] = (-dx)*dy*R5; tab[5] = (-dx)*dy*R5;
...@@ -361,15 +361,15 @@ public: ...@@ -361,15 +361,15 @@ public:
{ {
//Iteration over distantNeighbors //Iteration over distantNeighbors
int idxNeigh; int idxNeigh;
FPoint locCenter = local.getLeafCenter(); FPoint locCenter = getLeafCenter(local.getCoordinate());
for(idxNeigh=0 ; idxNeigh<343 ; idxNeigh++){ for(idxNeigh=0 ; idxNeigh<343 ; idxNeigh++){
//Need to test if current neighbor is one of the interaction list //Need to test if current neighbor is one of the interaction list
if(inIteractions[idxNeigh]){ if(distantNeighbors[idxNeigh]){
//Derivatives are computed iteratively //Derivatives are computed iteratively
FReal yetComputed[SizeVector]; //TODO C'est faux, k+n>p... FReal yetComputed[(2*P+1)*(2*P+2)*(2*P+3)/6];
FMemUtils::memset(yetComputed,0,SizeVector*sizeof(FReal(0))); FMemUtils::memset(yetComputed,0,SizeVector*sizeof(FReal(0)));
initDerivative(locCenter,inIteractions[idxNeigh].getLeafCenter(),yetComputed); initDerivative(locCenter,getLeafCenter(distantNeighbors[idxNeigh].getCoordinate),yetComputed);
//Iteration over Multipole / Local //Iteration over Multipole / Local
int al=0,bl=0,cl=0; int al=0,bl=0,cl=0;
...@@ -384,9 +384,9 @@ public: ...@@ -384,9 +384,9 @@ public:
//Iterating over multipole array : k //Iterating over multipole array : k
while(iterMult.hasNotFinished()){ while(iterMult.hasNotFinished()){
FReal data = computeDerivative(al+am,bl+bm,cl+cm,locCenter,inIteractions[idxNeigh].getLeafCenter(),yetComputed); FReal data = computeDerivative(al+am,bl+bm,cl+cm,locCenter,getLeafCenter(distantNeighbors[idxNeigh].getCoordinate()),yetComputed);
data *= iterMult.data()/fctl; data *= iterMult.data()/fctl;
iterLocal.setData(); iterLocal.setData(data);
//updating a,b,c and iterator //updating a,b,c and iterator
incPowers(&am,&bm,&cm); incPowers(&am,&bm,&cm);
...@@ -404,7 +404,7 @@ public: ...@@ -404,7 +404,7 @@ public:
/** /**
*@brief Divide and translate the local expansion of parent cell to child cell *@brief Translate the local expansion of parent cell to child cell
* *
*/ */
void L2L(const CellClass* const FRestrict local, void L2L(const CellClass* const FRestrict local,
...@@ -422,9 +422,43 @@ public: ...@@ -422,9 +422,43 @@ public:
void L2P(const CellClass* const local, void L2P(const CellClass* const local,
ContainerClass* const particles) ContainerClass* const particles)
{ {
} FPoint locCenter = getLeafCenter(local.getCoordinate());
//Iterator over particles
typename ContainerClass::ConstBasicIterator iterParticle(*particles);
while(iterParticle.hasNotFinished()){
FPoint dist = iterParticle.getPosition()-locCenter;
FReal dx = dist.getX();
FReal dy = dist.getY();
FReal dz = dist.getZ();
FReal forceX;
FReal forceY;
FReal forceZ;
int a=0,b=0,c=0;
//Iteration over Local array
typename ContainerClass::BasicIterator iterLocal(local.getLocal());
while(iterLocal.hasNotFinished()){
int idx = powerToIdx(a,b,c);
//Computation of forces
FReal locForce = iterLocal.data();
forceX = FMath::pow(dx,a)*locForce;
forceY = FMath::pow(dy,b)*locForce;
forceZ = FMath::pow(dz,c)*locForce;
//Application of forces
(iterParticle.data()).incForces(forceX,forceY,forceZ);
incPowers(&a,&b,&c);
iterLocal.gotoNext();
}
iterParticle.gotoNext();
}
}
}; };
#endif FTAYLORKERNEL_HPP #endif
...@@ -47,4 +47,4 @@ public: ...@@ -47,4 +47,4 @@ public:
}; };
#endif FTAYLORPARTICLE_HPP #endif
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