diff --git a/Src/Kernels/Taylor/FTaylorKernel.hpp b/Src/Kernels/Taylor/FTaylorKernel.hpp index 9e598929f2bcec56ae87234271f5328391e3d47a..d256bb70efc9c45a9104ec82949c160487549614 100644 --- a/Src/Kernels/Taylor/FTaylorKernel.hpp +++ b/Src/Kernels/Taylor/FTaylorKernel.hpp @@ -607,7 +607,7 @@ public: * where \f$x_c\f$ is the centre of the cell and \f$x_j\f$ the \f$j^{th}\f$ particles and \f$q_j\f$ its charge. */ void L2P(const CellClass* const local, - ContainerClass* const particles) + ContainerClass* const particles) { FPoint locCenter = getLeafCenter(local->getCoordinate()); @@ -631,13 +631,16 @@ public: //Iteration over particles for(int i=0 ; i<nbPart ; ++i){ - - FReal dx = posX[i]-locCenter.getX(); - FReal dy = posY[i]-locCenter.getY(); - FReal dz = posZ[i]-locCenter.getZ(); + // + FReal dx = posX[i] - locCenter.getX(); + FReal dy = posY[i] - locCenter.getY(); + FReal dz = posZ[i] - locCenter.getZ(); printf("dx = %f, dy = %f, dz = %f\n",dx,dy,dz); int a=0,b=0,c=0; + forceX[i] = 0.0 ; + forceY[i] = 0.0 ; + forceZ[i] = 0.0 ; // // Precompute an arrays of Array[i] = dx^(i-1) arrayDX[0] = 0.0 ; @@ -655,19 +658,24 @@ public: //Iteration over Local array const FReal * iterLocal = local->getLocal(); + FReal partPhyValue = phyValues[i]; + // + FReal locPot = 0.0, locForceX = 0.0, locForceY = 0.0, locForceZ = 0.0 ; for(int j=0 ; j<SizeVector ; ++j){ - - FReal partPhyValue = phyValues[i]; FReal locForce = iterLocal[j]; // compute the potential - targetsPotentials[i] += partPhyValue*iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1]; + locPot += iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1]; //Application of forces - forceX[i] -= partPhyValue*FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1]; - forceY[i] -= partPhyValue*FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1]; - forceZ[i] -= partPhyValue*FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c]; + locForceX += FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1]; + locForceY += FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1]; + locForceZ += FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c]; // incPowers(&a,&b,&c); } + targetsPotentials[i] = partPhyValue*locPot ; + forceX[i] -= partPhyValue*locForceX ; + forceY[i] -= partPhyValue*locForceY ; + forceZ[i] -= partPhyValue*locForceZ ; } } diff --git a/Tests/Kernels/testTaylorSimple.cpp b/Tests/Kernels/testTaylorSimple.cpp index 281f57e52d0e535149f7e563867c0b8672d59543..f57673d272ec2bbe84e67ef75c008323de3fdd2e 100644 --- a/Tests/Kernels/testTaylorSimple.cpp +++ b/Tests/Kernels/testTaylorSimple.cpp @@ -26,7 +26,7 @@ int main(int argc,char* argv[]){ static const int P = 3; static const int order = 1; - FPoint rootCenter(FReal(0),FReal(0),FReal(0)); + FPoint rootCenter(FReal(0.0),FReal(0.0),FReal(0.0)); FReal boxWidth = FReal(8); typedef FTaylorCell<P,order> CellClass; @@ -76,14 +76,21 @@ int main(int argc,char* argv[]){ for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ potential += potentials[idxPart]; - fx += forcesX[idxPart]; - fy += forcesY[idxPart]; - fz += forcesZ[idxPart]; + fx += forcesX[idxPart]; + fy += forcesY[idxPart]; + fz += forcesZ[idxPart]; printf("part : %f, fx = %f, fy = %f, fz = %f\n",leaf->getTargets()->getPositions()[0][0],fx,fy,fz); } }); - printf("potential : %f\n",potential); + // + FReal dx = part1Pos.getX() - part2Pos.getX(); + FReal dy = part1Pos.getY() - part2Pos.getY(); + FReal dz = part1Pos.getZ() - part2Pos.getZ(); + // + FReal potTheo = physVal2*physVal1*FMath::Sqrt(1.0 / (dx*dx + dy*dy + dz*dz)); + printf("Exact potential : %f Computed potential : %f Error: %f \n",potTheo, potential,potTheo- potential); + } return 0;