Commit 40f46768 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Change FComplexe to FComplex in all files and name of files

parent 077ce97b
......@@ -24,7 +24,7 @@
#include "../../Components/FBasicCell.hpp"
#include "../../Extensions/FExtendCellType.hpp"
#include "../../Utils/FComplexe.hpp"
#include "../../Utils/FComplex.hpp"
/**
* @author Pierre Blanchard (pierre.blanchard@inria.fr)
......@@ -51,17 +51,17 @@ class FInterpCell : public FBasicCell
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
// PB: Store multipole and local expansion in Fourier space
FComplexe transformed_multipole_exp[NRHS * NVALS * TransformedVectorSize];
FComplexe transformed_local_exp[NLHS * NVALS * TransformedVectorSize];
FComplex transformed_multipole_exp[NRHS * NVALS * TransformedVectorSize];
FComplex transformed_local_exp[NLHS * NVALS * TransformedVectorSize];
public:
FInterpCell(){
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
memset(transformed_multipole_exp, 0,
sizeof(FComplexe) * NRHS * NVALS * TransformedVectorSize);
sizeof(FComplex) * NRHS * NVALS * TransformedVectorSize);
memset(transformed_local_exp, 0,
sizeof(FComplexe) * NLHS * NVALS * TransformedVectorSize);
sizeof(FComplex) * NLHS * NVALS * TransformedVectorSize);
}
~FInterpCell() {}
......@@ -90,20 +90,20 @@ public:
}
/** Get Transformed Multipole */
const FComplexe* getTransformedMultipole(const int inRhs) const{
const FComplex* getTransformedMultipole(const int inRhs) const{
return this->transformed_multipole_exp + inRhs*TransformedVectorSize;
}
/** Get Transformed Local */
const FComplexe* getTransformedLocal(const int inRhs) const{
const FComplex* getTransformedLocal(const int inRhs) const{
return this->transformed_local_exp + inRhs*TransformedVectorSize;
}
/** Get Transformed Multipole */
FComplexe* getTransformedMultipole(const int inRhs){
FComplex* getTransformedMultipole(const int inRhs){
return this->transformed_multipole_exp + inRhs*TransformedVectorSize;
}
/** Get Transformed Local */
FComplexe* getTransformedLocal(const int inRhs){
FComplex* getTransformedLocal(const int inRhs){
return this->transformed_local_exp + inRhs*TransformedVectorSize;
}
......@@ -117,9 +117,9 @@ public:
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
memset(transformed_multipole_exp, 0,
sizeof(FComplexe) * NRHS * NVALS * TransformedVectorSize);
sizeof(FComplex) * NRHS * NVALS * TransformedVectorSize);
memset(transformed_local_exp, 0,
sizeof(FComplexe) * NLHS * NVALS * TransformedVectorSize);
sizeof(FComplex) * NLHS * NVALS * TransformedVectorSize);
}
///////////////////////////////////////////////////////
......@@ -171,7 +171,7 @@ public:
}
static constexpr int GetSize(){
return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplexe);
return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplex);
}
......
......@@ -34,10 +34,10 @@ namespace FP2P {
*/
template <typename MatrixKernelClass>
inline void MutualParticles(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue,
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
// Compute kernel of interaction...
const FPoint sourcePoint(sourceX,sourceY,sourceZ);
......@@ -75,9 +75,9 @@ inline void MutualParticles(const FReal sourceX,const FReal sourceY,const FReal
*/
template <typename MatrixKernelClass>
inline void NonMutualParticles(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
// Compute kernel of interaction...
const FPoint sourcePoint(sourceX,sourceY,sourceZ);
......@@ -112,7 +112,7 @@ inline void NonMutualParticles(const FReal sourceX,const FReal sourceY,const FRe
*/
template <class ContainerClass, typename MatrixKernelClass>
inline void FullMutualRIJ(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
const double CoreWidth2 = MatrixKernel->getCoreWidth2(); //PB: TODO directly call evaluateBlock
......@@ -122,158 +122,158 @@ inline void FullMutualRIJ(ContainerClass* const FRestrict inTargets, ContainerCl
const FReal*const targetsZ = inTargets->getPositions()[2];
for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials(i);
FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX(i);
FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY(i);
FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
sourcesForcesX[idxSource] -= coef[0];
sourcesForcesY[idxSource] -= coef[1];
sourcesForcesZ[idxSource] -= coef[2];
sourcesPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
}
}
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials(i);
FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX(i);
FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY(i);
FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
sourcesForcesX[idxSource] -= coef[0];
sourcesForcesY[idxSource] -= coef[1];
sourcesForcesZ[idxSource] -= coef[2];
sourcesPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
}
}
}
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
FReal dx = targetsX[idxSource] - targetsX[idxTarget];
FReal dy = targetsY[idxSource] - targetsY[idxTarget];
FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * targetsPhysicalValues[idxSource] );
targetsForcesX[idxSource] -= coef[0];
targetsForcesY[idxSource] -= coef[1];
targetsForcesZ[idxSource] -= coef[2];
targetsPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
for(int idxSource = idxTarget + 1 ; idxSource < nbParticlesTargets ; ++idxSource){
FReal dx = targetsX[idxSource] - targetsX[idxTarget];
FReal dy = targetsY[idxSource] - targetsY[idxTarget];
FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * targetsPhysicalValues[idxSource] );
targetsForcesX[idxSource] -= coef[0];
targetsForcesY[idxSource] -= coef[1];
targetsForcesZ[idxSource] -= coef[2];
targetsPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}// j
}// i
}
}
}
......@@ -282,7 +282,7 @@ inline void FullMutualRIJ(ContainerClass* const FRestrict inTargets, ContainerCl
*/
template <class ContainerClass, typename MatrixKernelClass>
inline void FullRemoteRIJ(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
const double CoreWidth2 = MatrixKernel->getCoreWidth2(); //PB: TODO directly call evaluateBlock
......@@ -292,79 +292,79 @@ inline void FullRemoteRIJ(ContainerClass* const FRestrict inTargets, ContainerCl
const FReal*const targetsZ = inTargets->getPositions()[2];
for(int idxNeighbors = 0 ; idxNeighbors < limiteNeighbors ; ++idxNeighbors){
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
// potential
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
}// j
}// i
}
}
}
for(int idxTarget = 0 ; idxTarget < nbParticlesTargets ; ++idxTarget){
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
// potential
FReal dx = sourcesX[idxSource] - targetsX[idxTarget];
FReal dy = sourcesY[idxSource] - targetsY[idxTarget];
FReal dz = sourcesZ[idxSource] - targetsZ[idxTarget];
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz + CoreWidth2);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal*const targetsPotentials = inTargets->getPotentials(i);
FReal*const targetsForcesX = inTargets->getForcesX(i);
FReal*const targetsForcesY = inTargets->getForcesY(i);
FReal*const targetsForcesZ = inTargets->getForcesZ(i);
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues(j);
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues(j);
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetsForcesX[idxTarget] += coef[0];
targetsForcesY[idxTarget] += coef[1];
targetsForcesZ[idxTarget] += coef[2];
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
}// j
}// i
}
}
}
}
}
......@@ -390,10 +390,10 @@ inline void FullRemoteRIJ(ContainerClass* const FRestrict inTargets, ContainerCl
*/
template<typename MatrixKernelClass>
inline void MutualParticlesRIJ(const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal* sourcePhysicalValue,
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal* targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential,
const FReal targetX,const FReal targetY,const FReal targetZ, const FReal* targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const MatrixKernelClass *const MatrixKernel){
const double CoreWidth2 = MatrixKernel->getCoreWidth2(); //PB: TODO directly call evaluateBlock
......@@ -409,52 +409,52 @@ inline void MutualParticlesRIJ(const FReal sourceX,const FReal sourceY,const FRe
FReal r[3]={dx,dy,dz};
for(unsigned int i = 0 ; i < 3 ; ++i){
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k)
coef[k]= -(targetPhysicalValue[j] * sourcePhysicalValue[j]);
// Grad of RIJ kernel is RIJK kernel => use same expression as in FInterpMatrixKernel
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
coef[k] *= FReal(3.) * ( FReal(-1.) + ri2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i=j!=k
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[k] * inv_distance_pow3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
coef[k] *= ( FReal(-1.) + FReal(3.) * ri2 * inv_distance_pow2 ) * r[j] * inv_distance_pow3;
else if(j==k) //i!=k=j
coef[k] *= ( FReal(-1.) + FReal(3.) * rj2 * inv_distance_pow2 ) * r[i] * inv_distance_pow3;
else //i!=k!=j
coef[k] *= FReal(3.) * r[i] * r[j] * r[k] * inv_distance_pow2 * inv_distance_pow3;
}
}// k
targetForceX[i] += coef[0];
targetForceY[i] += coef[1];
targetForceZ[i] += coef[2];
targetPotential[i] += ( potentialCoef * sourcePhysicalValue[j] );
sourceForceX[i] -= coef[0];
sourceForceY[i] -= coef[1];
sourceForceZ[i] -= coef[2];
sourcePotential[i] += ( potentialCoef * targetPhysicalValue[j] );
}// j
FReal ri2=r[i]*r[i];
for(unsigned int j = 0 ; j < 3 ; ++j){
// potentials
FReal potentialCoef;
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
potentialCoef = - r[i] * r[j] * inv_distance_pow3;
// forces
FReal rj2=r[j]*r[j];