Commit f814a7cf authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Provided generic functions for the evaluation of the direct interactions (pass...

Provided generic functions for the evaluation of the direct interactions (pass matrix kernel as argument).
parent 476eb362
......@@ -29,6 +29,7 @@
#include "Files/FFmaGenericLoader.hpp"
#include "Kernels/P2P/FP2P.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
//
......@@ -118,6 +119,9 @@ int main(int argc, char ** argv){
// ----------------------------------------------------------------------------------------------------------
// COMPUTATION
// ----------------------------------------------------------------------------------------------------------
// interaction kernel evaluator
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
FReal denergy = 0.0;
//
// computation
......@@ -142,7 +146,7 @@ int main(int argc, char ** argv){
particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
particles[idxTarget].getPosition().getZ(),particles[idxTarget].getPhysicalValue(),
&particles[idxTarget].setForces()[0],&particles[idxTarget].setForces()[1],
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential());
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),&MatrixKernel);
}
}
} // end for
......
......@@ -212,14 +212,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
}
void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/){
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
}
};
......
......@@ -462,14 +462,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles,AbstractBaseClass::MatrixKernel.getPtr());
}
void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/){
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27,AbstractBaseClass::MatrixKernel.getPtr());
}
};
......
......@@ -86,6 +86,7 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -94,11 +95,29 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
block[0] = one_over_r;
blockDerivative[0][0] = one_over_r3 * xy.getX();
blockDerivative[0][1] = one_over_r3 * xy.getY();
blockDerivative[0][2] = one_over_r3 * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -126,6 +145,7 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : LX(1.0),LY(1.0),LZ(1.0)
{ }
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -145,6 +165,23 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal one_over_rL = FReal(1.) / FMath::Sqrt(LX*xy.getX()*xy.getX() +
LY*xy.getY()*xy.getY() +
LZ*xy.getZ()*xy.getZ());
const FReal one_over_rL3 = one_over_rL*one_over_rL*one_over_rL;
block[0] = one_over_rL;
blockDerivative[0][0] = FMath::Sqrt(LX) * one_over_rL3 * xy.getX();
blockDerivative[0][1] = FMath::Sqrt(LY) * one_over_rL3 * xy.getY();
blockDerivative[0][2] = FMath::Sqrt(LZ) * one_over_rL3 * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -176,6 +213,7 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -184,11 +222,31 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
xy.getZ()*xy.getZ());
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal r2 = FReal(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
const FReal one_over_r2 = FReal(1.) / (r2);
const FReal one_over_r4 = one_over_r2*one_over_r2;
block[0] = one_over_r2;
const FReal coef = FReal(-2.) * one_over_r4;
blockDerivative[0][0] = coef * xy.getX();
blockDerivative[0][1] = coef * xy.getY();
blockDerivative[0][2] = coef * xy.getZ();
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -220,6 +278,7 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
int getPosition(const unsigned int) const
{return 0;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -231,11 +290,31 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
return one_over_r6 * one_over_r6 - one_over_r6;
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
// evaluate interaction and derivative (blockwise)
void evaluateBlockAndDerivative(const FPoint& x, const FPoint& y, FReal block[1], FReal blockDerivative[1][3]) const
{
const FPoint xy(x-y);
const FReal r = xy.norm();
const FReal r2 = r*r;
const FReal r3 = r2*r;
const FReal one_over_r6 = FReal(1.) / (r3*r3);
const FReal one_over_r8 = one_over_r6 / (r2);
block[0] = one_over_r6 * one_over_r6 - one_over_r6;
const FReal coef = FReal(12.0)*one_over_r6*one_over_r8 - FReal(6.0)*one_over_r8;
blockDerivative[0][0]= coef * xy.getX();
blockDerivative[0][1]= coef * xy.getY();
blockDerivative[0][2]= coef * xy.getZ();
}
FReal getScaleFactor(const FReal, const int) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
......@@ -317,6 +396,7 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
FReal getCoreWidth2() const
{return _CoreWidth2;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -343,6 +423,7 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
const FPoint xy(x-y);
......@@ -414,6 +495,7 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
FReal getCoreWidth2() const
{return _CoreWidth2;}
// evaluate interaction
FReal evaluate(const FPoint& x, const FPoint& y) const
{
// Convention for anti-symmetric kernels xy=y-x instead of xy=x-y
......@@ -460,6 +542,7 @@ struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
}
// evaluate interaction (blockwise)
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
// Convention for anti-symmetric kernels xy=y-x instead of xy=x-y
......
......@@ -4,69 +4,51 @@
#include "../P2P/FP2P.hpp"
template <KERNEL_FUNCTION_IDENTIFIER Identifier, int NVALS>
struct DirectInteractionComputer;
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteractionComputer<ONE_OVER_R, 1>
template <KERNEL_FUNCTION_IDENTIFIER Identifier, int NVALS>
struct DirectInteractionComputer
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const MatrixKernel){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14,MatrixKernel);
}
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
const int inSize,
const MatrixKernelClass *const MatrixKernel){
FP2P::FullRemote(inTargets,inNeighbors,inSize,MatrixKernel);
}
};
//! Specialization for Laplace potential on Non uniform domain
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteractionComputer<ONE_OVER_RH, 1>
struct DirectInteractionComputer<ONE_OVER_R, 1>
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const /*MatrixKernel*/){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
//
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
const int inSize,
const MatrixKernelClass *const /*MatrixKernel*/){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for Lennard-Jones potential */
template <>
struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, 1>
{
template <typename ContainerClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote(ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for GradGradR potential */
template <>
struct DirectInteractionComputer<R_IJ, 1>
......@@ -116,18 +98,20 @@ struct DirectInteractionComputer<R_IJK, 1>
template <int NVALS>
struct DirectInteractionComputer<ONE_OVER_R, NVALS>
{
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
ContainerClass* const NeighborSourceParticles[27],
const MatrixKernelClass *const /*MatrixKernel*/){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
template <typename ContainerClass, typename MatrixKernelClass>
static void P2PRemote(ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
const int inSize,
const MatrixKernelClass *const /*MatrixKernel*/){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
......@@ -135,29 +119,6 @@ struct DirectInteractionComputer<ONE_OVER_R, NVALS>
};
/*! Specialization for Lennard-Jones potential */
template <int NVALS>
struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, NVALS>
{
template <typename ContainerClass>
static void P2P(ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote(ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for GradGradR potential */
template <int NVALS>
struct DirectInteractionComputer<R_IJ, NVALS>
......
......@@ -10,84 +10,105 @@
*/
namespace FP2P {
/** P2P mutual interaction,
/**
* @brief MutualParticles (generic version)
* P2P mutual interaction,
* this function computes the interaction for 2 particles.
*
* Formulas are:
* \f[
* F = q_1 * q_2 / r^2
* P_1 = q_2 / r ; P_2 = q_1 / r
* F = - q_1 * q_2 * grad K{12}
* P_1 = q_2 * K{12} ; P_2 = q_1 * K_{12}
* \f]
* In details :
* \f$ F(x) = \frac{ \Delta_x * q_1 * q_2 }{ r^2 } = \Delta_x * F \f$
* In details for \f$K(x,y)=1/|x-y|=1/r\f$ :
* \f$ P_1 = \frac{ q_2 }{ r } \f$
* \f$ P_2 = \frac{ q_1 }{ r } \f$
* \f$ F(x) = \frac{ \Delta_x * q_1 * q_2 }{ r^2 } \f$
*
* @param sourceX
* @param sourceY
* @param sourceZ
* @param sourcePhysicalValue
* @param targetX
* @param targetY
* @param targetZ
* @param targetPhysicalValue
* @param targetForceX
* @param targetForceY
* @param targetForceZ
* @param targetPotential
* @param MatrixKernel pointer to an interaction kernel evaluator
*/
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
){
FReal dx = sourceX - targetX;
FReal dy = sourceY - targetY;
FReal dz = sourceZ - targetZ;
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_square_distance *= inv_distance;
inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
*targetForceX += dx;
*targetForceY += dy;
*targetForceZ += dz;
*targetPotential += ( inv_distance * sourcePhysicalValue );
*sourceForceX -= dx;
*sourceForceY -= dy;
*sourceForceZ -= dz;
*sourcePotential += ( inv_distance * targetPhysicalValue );
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);
const FPoint targetPoint(targetX,targetY,targetZ);
FReal Kxy[1];
FReal dKxy[1][3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint,targetPoint,Kxy,dKxy);
FReal coef = (targetPhysicalValue * sourcePhysicalValue);
(*targetForceX) += dKxy[0][0] * coef;
(*targetForceY) += dKxy[0][1] * coef;
(*targetForceZ) += dKxy[0][2] * coef;
(*targetPotential) += ( Kxy[0] * sourcePhysicalValue );
(*sourceForceX) -= dKxy[0][0] * coef;
(*sourceForceY) -= dKxy[0][1] * coef;
(*sourceForceZ) -= dKxy[0][2] * coef;
(*sourcePotential) += ( Kxy[0] * targetPhysicalValue );
}
/** P2P mutual interaction,
/**
* @brief NonMutualParticles (generic version)
* P2P mutual interaction,
* this function computes the interaction for 2 particles.
*
* Formulas are:
* \f[
* F = q_1 * q_2 / r^2
* P_1 = q_2 / r ; P_2 = q_1 / r
* F = - q_1 * q_2 * grad K{12}
* P_1 = q_2 * K{12} ; P_2 = q_1 * K_{12}
* \f]
* In details :
* \f$ F(x) = \frac{ \Delta_x * q_1 * q_2 }{ r^2 } = \Delta_x * F \f$
* In details for \f$K(x,y)=1/|x-y|=1/r\f$ :
* \f$ P_1 = \frac{ q_2 }{ r } \f$
* \f$ P_2 = \frac{ q_1 }{ r } \f$
* \f$ F(x) = \frac{ \Delta_x * q_1 * q_2 }{ r^2 } \f$
*/
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){
FReal dx = sourceX - targetX;
FReal dy = sourceY - targetY;
FReal dz = sourceZ - targetZ;
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
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);
const FPoint targetPoint(targetX,targetY,targetZ);
FReal Kxy[1];
FReal dKxy[1][3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint,targetPoint,Kxy,dKxy);
FReal coef = (targetPhysicalValue * sourcePhysicalValue);
(*targetForceX) += dKxy[0][0] * coef;
(*targetForceY) += dKxy[0][1] * coef;
(*targetForceZ) += dKxy[0][2] * coef;
(*targetPotential) += ( Kxy[0] * sourcePhysicalValue );
}
inv_square_distance *= inv_distance;
inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
*targetForceX += dx;
*targetForceY += dy;
*targetForceZ += dz;
*targetPotential += ( inv_distance * sourcePhysicalValue );
}
/*
* FullMutual (generic version)
*/
template <class ContainerClass, typename MatrixKernelClass>
inline void FullMutual(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors, const MatrixKernelClass *const MatrixKernel){
template <class ContainerClass>
inline void FullMutualLJ(ContainerClass* const FRestrict inTargets, ContainerClass* const inNeighbors[],
const int limiteNeighbors){
const int nbParticlesTargets = inTargets->getNbParticles();
const FReal*const targetsPhysicalValues = inTargets->getPhysicalValues();
......@@ -101,89 +122,72 @@ namespace FP2P {
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 sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX();
FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY();
FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ();
FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials();
if( inNeighbors[idxNeighbors] ){
const int nbParticlesSources = inNeighbors[idxNeighbors]->getNbParticles();
const FReal*const sourcesPhysicalValues = inNeighbors[idxNeighbors]->getPhysicalValues();
const FReal*const sourcesX = inNeighbors[idxNeighbors]->getPositions()[0];
const FReal*const sourcesY = inNeighbors[idxNeighbors]->getPositions()[1];
const FReal*const sourcesZ = inNeighbors[idxNeighbors]->getPositions()[2];
FReal*const sourcesForcesX = inNeighbors[idxNeighbors]->getForcesX();
FReal*const sourcesForcesY = inNeighbors[idxNeighbors]->getForcesY();
FReal*const sourcesForcesZ = inNeighbors[idxNeighbors]->getForcesZ();
FReal*const sourcesPotentials = inNeighbors[idxNeighbors]->getPotentials();
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];
for(int idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
FReal inv_distance_pow2 = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal inv_distance_pow6 = inv_distance_pow3 * inv_distance_pow3;
FReal inv_distance_pow8 = inv_distance_pow6 * inv_distance_pow2;
FReal coef = ((targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource])
* (FReal(12.0)*inv_distance_pow6*inv_distance_pow8 - FReal(6.0)*inv_distance_pow8));
FReal potentialCoef = (inv_distance_pow6*inv_distance_pow6-inv_distance_pow6);
dx *= coef;
dy *= coef;
dz *= coef;
targetsForcesX[idxTarget] += dx;
targetsForcesY[idxTarget] += dy;
targetsForcesZ[idxTarget] += dz;
targetsPotentials[idxTarget] += ( potentialCoef * sourcesPhysicalValues[idxSource] );
sourcesForcesX[idxSource] -= dx;
sourcesForcesY[idxSource] -= dy;
sourcesForcesZ[idxSource] -= dz;
sourcesPotentials[idxSource] += potentialCoef * targetsPhysicalValues[idxTarget];
}
}
// Compute kernel of interaction and its derivative
const FPoint sourcePoint(sourcesX[idxSource],sourcesY[idxSource],sourcesZ[idxSource]);
const FPoint targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
FReal Kxy[1];
FReal dKxy[1][3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint,targetPoint,Kxy,dKxy);
FReal coef = (targetsPhysicalValues[idxTarget] * sourcesPhysicalValues[idxSource]);
targetsForcesX[idxTarget] += dKxy[0][0] * coef;
targetsForcesY[idxTarget] += dKxy[0][1] * coef;
targetsForcesZ[idxTarget] += dKxy[0][2] * coef;
targetsPotentials[idxTarget] += ( Kxy[0] * sourcesPhysicalValues[idxSource] );
sourcesForcesX[idxSource] -= dKxy[0][0] * coef;
sourcesForcesY[idxSource] -= dKxy[0][1] * coef;
sourcesForcesZ[idxSource] -= dKxy[0][2] * coef;
sourcesPotentials[idxSource] += ( Kxy[0] * targetsPhysicalValues[idxTarget] );
}
}
}
}
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);
FReal inv_distance = FMath::Sqrt(inv_distance_pow2);
FReal inv_distance_pow3 = inv_distance_pow2 * inv_distance;
FReal inv_distance_pow6 = inv_distance_pow3 * inv_distance_pow3;
FReal inv_distance_pow8 = inv_distance_pow6 * inv_distance_pow2;
FReal coef = ((targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource])
* (FReal(12.0)*inv_distance_pow6*inv_distance_pow8 - FReal(6.0)*inv_distance_pow8));
FReal potentialCoef = (inv_distance_pow6*inv_distance_pow6-inv_distance_pow6);
dx *= coef;
dy *= coef;
dz *= coef;
targetsForcesX[idxTarget] += dx;