Commit c82fe260 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Major modification! Defined rigoreous target-source convention: p(x)=sum_y...

Major modification! Defined rigoreous target-source convention: p(x)=sum_y K(x,y)w(y) with x=target and y=source. Direct computation is done by calling MutualParticles(Target params,Source params) or P2P(Target, Source); Kernel matrices evaluation is done by calling evaluate(x=target,y=source) which returns k(x,y)=k(x-y).
parent 14a6b1c0
......@@ -134,13 +134,13 @@ int main(int argc, char ** argv){
//
for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
if( idxOther != idxTarget ){
FP2P::NonMutualParticles(
particles[idxOther].getPosition().getX(), particles[idxOther].getPosition().getY(),
particles[idxOther].getPosition().getZ(),particles[idxOther].getPhysicalValue(),
particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
FP2P::NonMutualParticles(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(),&MatrixKernel);
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),
particles[idxOther].getPosition().getX(), particles[idxOther].getPosition().getY(),
particles[idxOther].getPosition().getZ(),particles[idxOther].getPhysicalValue(),
&MatrixKernel);
}
}
} // end for
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FINTERPMATRIXKERNEL_TENSORIAL_INTERACTIONS_HPP
#define FINTERPMATRIXKERNEL_TENSORIAL_INTERACTIONS_HPP
#include <stdexcept>
#include "Utils/FPoint.hpp"
#include "Utils/FNoCopyable.hpp"
#include "Utils/FMath.hpp"
#include "Utils/FGlobal.hpp"
#include "FInterpMatrixKernel.hpp"
/**
*
* @author Pierre Blanchard (pierre.blanchard@inria.fr)
* @class FInterpMatrixKernels_TensorialInteractions
* Please read the license
*
* This class provides the evaluators and scaling functions of the matrix
* kernels. A matrix kernel should be understood in the sense of a kernel
* of interaction (or the fundamental solution of a given equation).
* It can either be scalar (NCMP=1) or tensorial (NCMP>1) depending on the
* dimension of the equation considered. NCMP denotes the number of components
* that are actually stored (e.g. 6 for a \f$3\times3\f$ symmetric tensor).
*
* Notes on the application scheme:
* Let there be a kernel \f$K\f$ such that \f$P_i(X)X=K_{ij}(X,Y)W_j(Y)\f$
* with \f$P\f$ the lhs of size NLHS and \f$W\f$ the rhs of size NRHS.
* The table applyTab provides the indices in the reduced storage table
* corresponding to the application scheme depicted earlier.
*
* PB: BEWARE! Homogeneous matrix kernels do not support cell width extension
* yet. Is it possible to find a reference width and a scale factor such that
* only 1 set of M2L opt can be used for all levels??
*
*/
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//
// Tensorial Matrix Kernels (NCMP>1)
//
// The definition of the potential p and force f are extended to the case
// of tensorial interaction kernels:
// p_i(x) = K_{ip}(x,y)w_p(y), \forall i=1..NPOT, p=1..NPV
// f_{ik}= w_p(x)K_{ip,k}(x,y)w_p(y) "
//
// Since the interpolation scheme is such that
// p_i(x) \approx S^m(x) L^{m}_{ip}
// f_{ik}= w_p(x) \nabla_k S^m(x) L^{m}_{ip}
// with
// L^{m}_{ip} = K^{mn}_{ip} S^n(y) w_p(y) (local expansion)
// M^{m}_{p} = S^n(y) w_p(y) (multipole expansion)
// then the multipole exp have NPV components and the local exp NPOT*NPV.
//
// NB1: Only the computation of forces requires that the sum over p is
// performed at L2P step. It could be done at M2L step for the potential.
//
// NB2: An efficient application of the matrix kernel is highly kernel
// dependent, we recommand overriding the P2M/M2L/L2P function of the kernel
// you are using in order to have opsimal performances + set your own NRHS/NLHS.
//
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/// R_{,ij}
// PB: IMPORTANT! This matrix kernel does not present the symmetries
// required by ChebSym kernel => only suited for Unif interpolation
template <class FReal>
struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel<FReal>
{
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const unsigned int NK = 3*3; //< total number of components
static const unsigned int NCMP = 6; //< number of independent components
static const unsigned int NPV = 3; //< dim of physical values
static const unsigned int NPOT = 3; //< dim of potentials
static const unsigned int NRHS = NPV; //< dim of mult exp
static const unsigned int NLHS = NPOT*NPV; //< dim of loc exp (no csummation over j during M2L)
// store indices (i,j) corresponding to sym idx
static const unsigned int indexTab[/*2*NCMP=12*/];
// store positions in sym tensor (when looping over NRHSxNLHS)
static const unsigned int applyTab[/*NK=9*/];
// indices to be set at construction if component-wise evaluation is performed
const unsigned int _i,_j;
// Material Parameters
const FReal _CoreWidth2; // if >0 then kernel is NON homogeneous
FInterpMatrixKernel_R_IJ(const FReal CoreWidth = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+NCMP]), _CoreWidth2(CoreWidth*CoreWidth)
{}
// copy ctor
FInterpMatrixKernel_R_IJ(const FInterpMatrixKernel_R_IJ& other)
: _i(other._i), _j(other._j), _CoreWidth2(other._CoreWidth2)
{}
static const char* getID() { return "R_IJ"; }
static void printInfo() {
std::cout << "K_ij(x,y)=r_{,ij}=d^2/dx_idx_j(r), with r=|x-y|_a=sqrt(a^2 + (x_i-y_i)^2)." << std::endl;
std::cout << "Compute potentials p_i(x)=sum_{y,j}(K_ij(x,y)w_j(y))" << std::endl;
std::cout << " forces f_{ik}(x)=-d/dx_k(sum_{y,j}(w_j(x)K_ij(x,y)w_j(y)))" << std::endl;
std::cout << std::endl;
}
// returns position in reduced storage from position in full 3x3 matrix
unsigned int getPosition(const unsigned int n) const
{return applyTab[n];}
// returns Core Width squared
FReal getCoreWidth2() const
{return _CoreWidth2;}
// returns coefficient of mutual interaction
// 1 for symmetric kernels
// -1 for antisymmetric kernels
// somethings else if other property of symmetry
FReal getMutualCoefficient() const{ return FReal(1.); }
// evaluate interaction
template <class ValueClass>
FReal evaluate(const ValueClass& xt, const ValueClass& yt, const ValueClass& zt,
const ValueClass& xs, const ValueClass& ys, const ValueClass& zs) const
{
const ValueClass diffx = (xt-xs);
const ValueClass diffy = (yt-ys);
const ValueClass diffz = (zt-zs);
const ValueClass r2 = diffx*diffx+diffy*diffy+diffz*diffz;
const ValueClass one_over_r = FMath::One<ValueClass>()/FMath::Sqrt(r2 + FMath::ConvertTo<ValueClass,FReal>(_CoreWidth2));
const ValueClass one_over_r3 = one_over_r*one_over_r*one_over_r;
ValueClass ri,rj;
if(_i==0) ri=diffx;
else if(_i==1) ri=diffy;
else if(_i==2) ri=diffz;
else throw std::runtime_error("Update i!");
if(_j==0) rj=diffx;
else if(_j==1) rj=diffy;
else if(_j==2) rj=diffz;
else throw std::runtime_error("Update j!");
if(_i==_j)
return one_over_r - ri * ri * one_over_r3;
else
return - ri * rj * one_over_r3;
}
// evaluate interaction (blockwise)
template <class ValueClass>
void evaluateBlock(const ValueClass& xt, const ValueClass& yt, const ValueClass& zt,
const ValueClass& xs, const ValueClass& ys, const ValueClass& zs,
ValueClass* block) const
{
const ValueClass diffx = (xt-xs);
const ValueClass diffy = (yt-ys);
const ValueClass diffz = (zt-zs);
const ValueClass r2 = diffx*diffx+diffy*diffy+diffz*diffz;
const ValueClass one_over_r = FMath::One<ValueClass>()/FMath::Sqrt(r2 + FMath::ConvertTo<ValueClass,FReal>(_CoreWidth2));
const ValueClass one_over_r3 = one_over_r*one_over_r*one_over_r;
const ValueClass r[3] = {diffx,diffy,diffz};
for(unsigned int d=0;d<NCMP;++d){
unsigned int i = indexTab[d];
unsigned int j = indexTab[d+NCMP];
if(i==j)
block[d] = one_over_r - r[i] * r[i] * one_over_r3;
else
block[d] = - r[i] * r[j] * one_over_r3;
}
}
// evaluate interaction and derivative (blockwise)
// [TODO] Fix! Add corewidth!
template <class ValueClass>
void evaluateBlockAndDerivative(const ValueClass& xt, const ValueClass& yt, const ValueClass& zt,
const ValueClass& xs, const ValueClass& ys, const ValueClass& zs,
ValueClass block[NCMP], ValueClass blockDerivative[NCMP][3]) const
{
const ValueClass diffx = (xt-xs);
const ValueClass diffy = (yt-ys);
const ValueClass diffz = (zt-zs);
const ValueClass r2[3] = {diffx*diffx,diffy*diffy,diffz*diffz};
const ValueClass one_over_r2 = FMath::One<ValueClass>() / (r2[0] + r2[1] + r2[2] + FMath::ConvertTo<ValueClass,FReal>(_CoreWidth2));
const ValueClass one_over_r = FMath::Sqrt(one_over_r2);
const ValueClass one_over_r3 = one_over_r2*one_over_r;
const ValueClass r[3] = {diffx,diffy,diffz};
const ValueClass Three = FMath::ConvertTo<ValueClass,FReal>(3.);
const ValueClass MinusOne = - FMath::One<ValueClass>();
for(unsigned int d=0;d<NCMP;++d){
unsigned int i = indexTab[d];
unsigned int j = indexTab[d+NCMP];
// evaluate kernel
if(i==j)
block[d] = one_over_r - r2[i] * one_over_r3;
else
block[d] = - r[i] * r[j] * one_over_r3;
// evaluate derivative
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
if(j==k) //i=j=k
blockDerivative[d][k] = Three * ( MinusOne + r2[i] * one_over_r2 ) * r[i] * one_over_r3;
else //i=j!=k
blockDerivative[d][k] = ( MinusOne + Three * r2[i] * one_over_r2 ) * r[k] * one_over_r3;
}
else{ //(i!=j)
if(i==k) //i=k!=j
blockDerivative[d][k] = ( MinusOne + Three * r2[i] * one_over_r2 ) * r[j] * one_over_r3;
else if(j==k) //i!=k=j
blockDerivative[d][k] = ( MinusOne + Three * r2[j] * one_over_r2 ) * r[i] * one_over_r3;
else //i!=k!=j
blockDerivative[d][k] = Three * r[i] * r[j] * r[k] * one_over_r2 * one_over_r3;
}
}// k
}// NCMP
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
// R_{,ij} is homogeneous to [L]/[L]^{-2}=[L]^{-1}
// => scales like ONE_OVER_R
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(2.) / CellWidth;
}
FReal evaluate(const FPoint<FReal>& pt, const FPoint<FReal>& ps) const{
return evaluate<FReal>(pt.getX(), pt.getY(), pt.getZ(), ps.getX(), ps.getY(), ps.getZ());
}
void evaluateBlock(const FPoint<FReal>& pt, const FPoint<FReal>& ps, FReal* block) const{
evaluateBlock<FReal>(pt.getX(), pt.getY(), pt.getZ(), ps.getX(), ps.getY(), ps.getZ(), block);
}
void evaluateBlockAndDerivative(const FPoint<FReal>& pt, const FPoint<FReal>& ps,
FReal block[NCMP], FReal blockDerivative[NCMP][3]) const {
evaluateBlockAndDerivative<FReal>(pt.getX(), pt.getY(), pt.getZ(), ps.getX(), ps.getY(), ps.getZ(), block, blockDerivative);
}
};
/// R_IJ
template <class FReal>
const unsigned int FInterpMatrixKernel_R_IJ<FReal>::indexTab[]={0,0,0,1,1,2,
0,1,2,1,2,2};
template <class FReal>
const unsigned int FInterpMatrixKernel_R_IJ<FReal>::applyTab[]={0,1,2,
1,3,4,
2,4,5};
#endif // FINTERPMATRIXKERNEL_TENSORIAL_INTERACTIONS_HPP
// [--END--]
This diff is collapsed.
......@@ -58,8 +58,10 @@ namespace FP2P {
const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
FReal Kxy[1];
FReal dKxy[3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
MatrixKernel->evaluateBlockAndDerivative(targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),
sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
Kxy,dKxy);
const FReal mutual_coeff = MatrixKernel->getMutualCoefficient(); // 1 if symmetric; -1 if antisymmetric
for(int idxVals = 0 ; idxVals < NVALS ; ++idxVals){
......@@ -76,7 +78,7 @@ namespace FP2P {
sourcesForcesX[idxSourceValue] -= dKxy[0] * coef;
sourcesForcesY[idxSourceValue] -= dKxy[1] * coef;
sourcesForcesZ[idxSourceValue] -= dKxy[2] * coef;
sourcesPotentials[idxSourceValue] += ( Kxy[0] * targetsPhysicalValues[idxTargetValue] );
sourcesPotentials[idxSourceValue] += ( mutual_coeff * Kxy[0] * targetsPhysicalValues[idxTargetValue] );
} // NVALS
......@@ -109,8 +111,10 @@ namespace FP2P {
const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
FReal Kxy[1];
FReal dKxy[3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
MatrixKernel->evaluateBlockAndDerivative(targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),
sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
Kxy,dKxy);
const FReal mutual_coeff = MatrixKernel->getMutualCoefficient(); // 1 if symmetric; -1 if antisymmetric
for(int idxVals = 0 ; idxVals < NVALS ; ++idxVals){
......@@ -127,7 +131,7 @@ namespace FP2P {
targetsForcesX[idxSourceValue] -= dKxy[0] * coef;
targetsForcesY[idxSourceValue] -= dKxy[1] * coef;
targetsForcesZ[idxSourceValue] -= dKxy[2] * coef;
targetsPotentials[idxSourceValue] += ( Kxy[0] * targetsPhysicalValues[idxTargetValue] );
targetsPotentials[idxSourceValue] += ( mutual_coeff * Kxy[0] * targetsPhysicalValues[idxTargetValue] );
}// NVALS
......@@ -172,8 +176,9 @@ inline void FullRemoteMultiRhs(ContainerClass* const FRestrict inTargets, const
const FPoint<FReal> targetPoint(targetsX[idxTarget],targetsY[idxTarget],targetsZ[idxTarget]);
FReal Kxy[1];
FReal dKxy[3];
MatrixKernel->evaluateBlockAndDerivative(sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),Kxy,dKxy);
MatrixKernel->evaluateBlockAndDerivative(targetPoint.getX(),targetPoint.getY(),targetPoint.getZ(),
sourcePoint.getX(),sourcePoint.getY(),sourcePoint.getZ(),
Kxy,dKxy);
for(int idxVals = 0 ; idxVals < NVALS ; ++idxVals){
......
......@@ -25,14 +25,13 @@
*/
namespace FP2PR{
template <class FReal>
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;
inline void MutualParticles(const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue,
FReal* sourceForceX, FReal* sourceForceY, FReal* sourceForceZ, FReal* sourcePotential){
FReal dx = targetX - sourceX;
FReal dy = targetY - sourceY;
FReal dz = targetZ - sourceZ;
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
......@@ -40,9 +39,9 @@ inline void MutualParticles(const FReal sourceX,const FReal sourceY,const FReal
inv_square_distance *= inv_distance;
inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
*targetForceX += dx;
*targetForceY += dy;
......@@ -56,12 +55,12 @@ inline void MutualParticles(const FReal sourceX,const FReal sourceY,const FReal
}
template <class FReal>
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;
inline void NonMutualParticles(const FReal targetX,const FReal targetY,const FReal targetZ, const FReal targetPhysicalValue,
FReal* targetForceX, FReal* targetForceY, FReal* targetForceZ, FReal* targetPotential,
const FReal sourceX,const FReal sourceY,const FReal sourceZ, const FReal sourcePhysicalValue){
FReal dx = targetX - sourceX;
FReal dy = targetY - sourceY;
FReal dz = targetZ - sourceZ;
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
FReal inv_distance = FMath::Sqrt(inv_square_distance);
......@@ -69,9 +68,10 @@ inline void NonMutualParticles(const FReal sourceX,const FReal sourceY,const FRe
inv_square_distance *= inv_distance;
inv_square_distance *= targetPhysicalValue * sourcePhysicalValue;
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
// d/dx(1/|x-y|)=-(x-y)/r^3
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
*targetForceX += dx;
*targetForceY += dy;
......@@ -119,9 +119,9 @@ static void GenericFullMutual(ContainerClass* const FRestrict inTargets, Contain
ComputeClass tpo = FMath::Zero<ComputeClass>();
for(FSize idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
ComputeClass dx = sourcesX[idxSource] - tx;
ComputeClass dy = sourcesY[idxSource] - ty;
ComputeClass dz = sourcesZ[idxSource] - tz;
ComputeClass dx = tx - sourcesX[idxSource];
ComputeClass dy = ty - sourcesY[idxSource];
ComputeClass dz = tz - sourcesZ[idxSource];
ComputeClass inv_square_distance = mOne / (dx*dx + dy*dy + dz*dz);
const ComputeClass inv_distance = FMath::Sqrt(inv_square_distance);
......@@ -129,9 +129,9 @@ static void GenericFullMutual(ContainerClass* const FRestrict inTargets, Contain
inv_square_distance *= inv_distance;
inv_square_distance *= tv * sourcesPhysicalValues[idxSource];
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
tfx += dx;
tfy += dy;
......@@ -193,18 +193,18 @@ static void GenericInner(ContainerClass* const FRestrict inTargets){
for(FSize idxSource = (idxTarget+NbFRealInComputeClass)/NbFRealInComputeClass ; idxSource < nbParticlesSources ; ++idxSource){
ComputeClass dx = sourcesX[idxSource] - tx;
ComputeClass dy = sourcesY[idxSource] - ty;
ComputeClass dz = sourcesZ[idxSource] - tz;
ComputeClass dx = tx - sourcesX[idxSource];
ComputeClass dy = ty - sourcesY[idxSource];
ComputeClass dz = tz - sourcesZ[idxSource];
ComputeClass inv_square_distance = mOne / (dx*dx + dy*dy + dz*dz);
const ComputeClass inv_distance = FMath::Sqrt(inv_square_distance);
inv_square_distance *= inv_distance;
inv_square_distance *= tv * sourcesPhysicalValues[idxSource];
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
tfx += dx;
tfy += dy;
......@@ -228,9 +228,9 @@ static void GenericInner(ContainerClass* const FRestrict inTargets){
const FSize limitForTarget = NbFRealInComputeClass-(idxTarget%NbFRealInComputeClass);
for(FSize idxS = 1 ; idxS < limitForTarget ; ++idxS){
const FSize idxSource = idxTarget + idxS;
FReal dx = targetsX[idxSource] - targetsX[idxTarget];
FReal dy = targetsY[idxSource] - targetsY[idxTarget];
FReal dz = targetsZ[idxSource] - targetsZ[idxTarget];
FReal dx = targetsX[idxTarget] - targetsX[idxSource];
FReal dy = targetsY[idxTarget] - targetsY[idxSource];
FReal dz = targetsZ[idxTarget] - targetsZ[idxSource];
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz);
const FReal inv_distance = FMath::Sqrt(inv_square_distance);
......@@ -238,9 +238,9 @@ static void GenericInner(ContainerClass* const FRestrict inTargets){
inv_square_distance *= inv_distance;
inv_square_distance *= targetsPhysicalValues[idxTarget] * targetsPhysicalValues[idxSource];
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
targetsForcesX[idxTarget] += dx;
targetsForcesY[idxTarget] += dy;
......@@ -289,9 +289,9 @@ static void GenericFullRemote(ContainerClass* const FRestrict inTargets, const C
ComputeClass tpo = FMath::Zero<ComputeClass>();
for(FSize idxSource = 0 ; idxSource < nbParticlesSources ; ++idxSource){
ComputeClass dx = sourcesX[idxSource] - tx;
ComputeClass dy = sourcesY[idxSource] - ty;
ComputeClass dz = sourcesZ[idxSource] - tz;
ComputeClass dx = tx - sourcesX[idxSource];
ComputeClass dy = ty - sourcesY[idxSource];
ComputeClass dz = tz - sourcesZ[idxSource];
ComputeClass inv_square_distance = mOne / (dx*dx + dy*dy + dz*dz);
const ComputeClass inv_distance = FMath::Sqrt(inv_square_distance);
......@@ -299,9 +299,9 @@ static void GenericFullRemote(ContainerClass* const FRestrict inTargets, const C
inv_square_distance *= inv_distance;
inv_square_distance *= tv * sourcesPhysicalValues[idxSource];
dx *= inv_square_distance;
dy *= inv_square_distance;
dz *= inv_square_distance;
dx *= - inv_square_distance;
dy *= - inv_square_distance;
dz *= - inv_square_distance;
tfx += dx;
tfy += dy;
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FP2P_TENSORIALKIJ_HPP
#define FP2P_TENSORIALKIJ_HPP
namespace FP2P {
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
// Tensorial Matrix Kernels: K_IJ / p_i=\sum_j K_{ij} w_j
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
/**
* @brief MutualParticlesKIJ
* @param sourceX
* @param sourceY
* @param sourceZ
* @param sourcePhysicalValue
* @param sourceForceX
* @param sourceForceY
* @param sourceForceZ
* @param sourcePotential
* @param targetX
* @param targetY
* @param targetZ
* @param targetPhysicalValue
* @param targetForceX
* @param targetForceY
* @param targetForceZ
* @param targetPotential
*/
template<class FReal, typename MatrixKernelClass>
inline void MutualParticlesKIJ(const FReal targetX