FAbstractSphericalKernel.hpp 32.5 KB
Newer Older
1
// See LICENCE file at project root
2 3
#ifndef FABSTRACTSPHERICALKERNEL_HPP
#define FABSTRACTSPHERICALKERNEL_HPP
4

5
#include <iostream>
6
#include "../../Components/FAbstractKernels.hpp"
7

8
#include "../../Utils/FGlobal.hpp"
9

10 11
#include "../../Utils/FMemUtils.hpp"
#include "../../Utils/FSmartPointer.hpp"
12
#include "../../Utils/FPoint.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
13
#include "../../Utils/FAssert.hpp"
14

15
#include "../../Containers/FTreeCoordinate.hpp"
16

17
#include "../P2P/FP2PR.hpp"
18

19 20 21 22
#include "FHarmonic.hpp"

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
23
* This is the abstract spherical harmonic kernel
24
*/
25
template< class FReal, class CellClass, class ContainerClass>
26
class FAbstractSphericalKernel : public FAbstractKernels<CellClass,ContainerClass> {
27
protected:
28 29 30
    const int   devP;           //< The P
    const FReal boxWidth;       //< the box width at leaf level
    const int   treeHeight;     //< The height of the tree
31

32 33
    const FReal widthAtLeafLevel;       //< the width of a box at leaf level
    const FReal widthAtLeafLevelDiv2;   //< the width of a box at leaf level divided by 2
34
    const FPoint<FReal> boxCorner;        //< the corner of the box system
35

36
    FHarmonic<FReal> harmonic; //< The harmonic computation class
37

38
    // For normal computation
39 40
    FSmartPointer<FComplex<FReal>*> preL2LTransitions; //< The pre-computation for the L2L based on the level
    FSmartPointer<FComplex<FReal>*> preM2MTransitions; //< The pre-computation for the M2M based on the level
41

42 43

    /** Alloc and init pre-vectors*/
44
    void allocAndInit(){
45 46 47 48
        preL2LTransitions = new FComplex<FReal>*[treeHeight ];
        memset(preL2LTransitions.getPtr(), 0, (treeHeight) * sizeof(FComplex<FReal>*));
        preM2MTransitions = new FComplex<FReal>*[treeHeight];
        memset(preM2MTransitions.getPtr(), 0, (treeHeight) * sizeof(FComplex<FReal>*));
49

50 51
        FReal treeWidthAtLevel = (boxWidth)/2;
        for(int idxLevel = 0 ; idxLevel < treeHeight - 1 ; ++idxLevel ){
52 53
            preL2LTransitions[idxLevel] = new FComplex<FReal>[ 8 * harmonic.getExpSize()];
            preM2MTransitions[idxLevel] = new FComplex<FReal>[ 8 * harmonic.getExpSize()];
54

55
            const FPoint<FReal> father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel);
56 57 58 59
            treeWidthAtLevel /= 2;

            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild ){
                FTreeCoordinate childBox;
60
                childBox.setPositionFromMorton(idxChild);
61

62
                const FPoint<FReal> M2MVector (
63 64 65 66 67
                        father.getX() - (treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))),
                        father.getY() - (treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))),
                        father.getZ() - (treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2)))
                        );

68 69
                harmonic.computeInner(FSpherical<FReal>(M2MVector));
                FMemUtils::copyall<FComplex<FReal>>(&preM2MTransitions[idxLevel][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
70

71
                const FPoint<FReal> L2LVector (
72 73 74 75 76
                        (treeWidthAtLevel * FReal(1 + (childBox.getX() * 2))) - father.getX(),
                        (treeWidthAtLevel * FReal(1 + (childBox.getY() * 2))) - father.getY(),
                        (treeWidthAtLevel * FReal(1 + (childBox.getZ() * 2))) - father.getZ()
                        );

77 78
                harmonic.computeInner(FSpherical<FReal>(L2LVector));
                FMemUtils::copyall<FComplex<FReal>>(&preL2LTransitions[idxLevel][harmonic.getExpSize() * idxChild], harmonic.result(), harmonic.getExpSize());
79 80
           }
        }
81 82
    }

83
    /** Get a leaf real position from its tree coordinate */
84 85
    FPoint<FReal> getLeafCenter(const FTreeCoordinate coordinate) const {
        return FPoint<FReal>(
86 87 88 89 90
                    FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
                    FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
                    FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX());
    }

91

92
public:
93
    /** Kernel constructor */
94
    FAbstractSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter)
95 96 97
        : devP(inDevP),
          boxWidth(inBoxWidth),
          treeHeight(inTreeHeight),
98
          widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))),
99 100 101
          widthAtLeafLevelDiv2(widthAtLeafLevel/2),
          boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),inBoxCenter.getZ()-(inBoxWidth/2)),
          harmonic(inDevP),
102 103
          preL2LTransitions(nullptr),
          preM2MTransitions(nullptr) {
104

105 106 107 108
        allocAndInit();
    }

    /** Copy constructor */
109
    FAbstractSphericalKernel(const FAbstractSphericalKernel& other)
110 111 112 113 114 115 116
        : devP(other.devP),
          boxWidth(other.boxWidth),
          treeHeight(other.treeHeight),
          widthAtLeafLevel(other.widthAtLeafLevel),
          widthAtLeafLevelDiv2(other.widthAtLeafLevelDiv2),
          boxCorner(other.boxCorner),
          harmonic(other.devP),
117 118
          preL2LTransitions(other.preL2LTransitions),
          preM2MTransitions(other.preM2MTransitions) {
119

120 121 122
    }

    /** Default destructor */
123
    virtual ~FAbstractSphericalKernel(){
124
        if(preL2LTransitions.isLast()){
125
            FMemUtils::DeleteAllArray(preL2LTransitions.getPtr(), treeHeight);
126 127
        }
        if(preM2MTransitions.isLast()){
128
            FMemUtils::DeleteAllArray(preM2MTransitions.getPtr(), treeHeight);
129
        }
130 131 132
    }

    /** P2M with a cell and all its particles */
133 134 135 136 137 138
    template<class SymbolicData>
    void P2M(typename CellClass::multipole_t* const LeafMultipole,
             const SymbolicData* const LeafSymb,
             const ContainerClass* const inParticles)
    {
        FComplex<FReal>* FRestrict const cellMultiPole = LeafMultipole->get();
139
        // Copying the position is faster than using cell position
140
        const FPoint<FReal> polePosition = getLeafCenter(LeafSymb->getCoordinate());
141
        // For all particles in the leaf box
142 143 144 145
        const FReal*const physicalValues = inParticles->getPhysicalValues();
        const FReal*const positionsX = inParticles->getPositions()[0];
        const FReal*const positionsY = inParticles->getPositions()[1];
        const FReal*const positionsZ = inParticles->getPositions()[2];
146
        for(FSize idxPart = 0 ; idxPart < inParticles->getNbParticles() ; ++idxPart){
147
            // P2M
148
            particleToMultiPole(cellMultiPole, polePosition,
149
                                FPoint<FReal>(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]),
150
                                physicalValues[idxPart]);
151 152 153
        }
    }

154
    /** M2M with a cell and all its child */
155 156 157 158 159 160 161 162
    template<class SymbolicData>
    void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole,
             const SymbolicData* const ParentSymb,
             const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles,
             const SymbolicData* const /*ChildSymbs*/[])
    {
        int level = static_cast<int>(ParentSymb->getLevel());
        FComplex<FReal>* FRestrict const multipole_exp_target = ParentMultipole->get();
163
        // iter on each child and process M2M
164 165
        const FComplex<FReal>* FRestrict const preM2MTransitionsAtLevel
            = preM2MTransitions[level];
166
        for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
167 168
            if(ChildMultipoles[idxChild]){
                multipoleToMultipole(multipole_exp_target, ChildMultipoles[idxChild]->get(), &preM2MTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
169 170
            }
        }
171 172
    }

173
    /** M2L with a cell and all the existing neighbors */
174 175 176
    // virtual void M2L(CellClass* const FRestrict inLocal, const CellClass* inInteractions[],
    //                  const int neighborPositions[], const int inSize, const int inLevel) override
    //     = 0;
177

178
    /** L2L with a cell and all its child */
179 180 181 182 183 184 185 186 187 188
    template<class SymbolicData>
    void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion,
             const SymbolicData* const ParentSymb,
             typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansions,
             const SymbolicData* const /*ChildSymbs*/[])
    {
        int level = static_cast<int>(ParentSymb->getLevel());
        // iter on each ChildrenExpansions and process L2L
        const FComplex<FReal>* FRestrict const preL2LTransitionsAtLevel
            = preL2LTransitions[level];
189
        for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
190 191 192 193
            if(ChildExpansions[idxChild]){
                localToLocal(ChildExpansions[idxChild]->get(),
                             ParentExpansion->get(),
                             &preL2LTransitionsAtLevel[idxChild * harmonic.getExpSize()]);
194 195
            }
        }
196 197 198
    }

    /** L2P with a cell and all its particles */
199 200 201 202 203 204
    template<class SymbolicData>
    void L2P(const typename CellClass::local_expansion_t* const local,
             const SymbolicData* const LeafSymb,
             ContainerClass* const inParticles)
    {
        const FComplex<FReal>* const cellLocal = local->get();
205
        // Copying the position is faster than using cell position
206
        const FPoint<FReal> localPosition = getLeafCenter(LeafSymb->getCoordinate());
207
        // For all particles in the leaf box
208 209 210 211 212 213 214 215
        const FReal*const physicalValues = inParticles->getPhysicalValues();
        const FReal*const positionsX = inParticles->getPositions()[0];
        const FReal*const positionsY = inParticles->getPositions()[1];
        const FReal*const positionsZ = inParticles->getPositions()[2];
        FReal*const potentials = inParticles->getPotentials();
        FReal*const forcesX = inParticles->getForcesX();
        FReal*const forcesY = inParticles->getForcesY();
        FReal*const forcesZ = inParticles->getForcesZ();
216
        for(FSize idxPart = 0 ; idxPart < inParticles->getNbParticles() ; ++idxPart){
217
            // L2P
218
            localToParticle(localPosition, cellLocal,
219
                            FPoint<FReal>(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]),
220 221
                            physicalValues[idxPart], &potentials[idxPart],
                            &forcesX[idxPart],&forcesY[idxPart],&forcesZ[idxPart]);
222 223 224
        }
    }

225 226 227 228 229
    /** P2P
      * This function proceed the P2P using particlesMutualInteraction
      * The computation is done for interactions with an index <= 13.
      * (13 means current leaf (x;y;z) = (0;0;0)).
      * Calling this method in multi thread should be done carrefully.
230
      */
231
    void P2P(const FTreeCoordinate& inPosition,
232
             ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict inSources,
233 234
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
235 236 237 238 239 240 241 242 243
        if(inTargets == inSources){
            FP2PRT<FReal>::template Inner<ContainerClass>(inTargets);
            P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
        }
        else{
            const ContainerClass* const srcPtr[1] = {inSources};
            FP2PRT<FReal>::template FullRemote<ContainerClass>(inTargets,srcPtr,1);
            FP2PRT<FReal>::template FullRemote<ContainerClass>(inTargets,inNeighbors,inSize);
        }
244 245 246 247 248 249
    }

    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
             ContainerClass* const FRestrict inTargets,
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
250 251 252 253 254 255
        int nbNeighborsToCompute = 0;
        while(nbNeighborsToCompute < inSize
              && neighborPositions[nbNeighborsToCompute] < 14){
            nbNeighborsToCompute += 1;
        }
        FP2PRT<FReal>::template FullMutual<ContainerClass>(inTargets,inNeighbors,nbNeighborsToCompute);
berenger-bramas's avatar
berenger-bramas committed
256 257
    }

258 259 260 261

    /** Use mutual even if it not useful and call particlesMutualInteraction */
    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
262
                   const ContainerClass* const inNeighbors[], const int neighborPositions[],
263 264
                   const int inSize) override {
        FP2PRT<FReal>::template FullRemote<ContainerClass>(inTargets,inNeighbors,inSize);
265 266
    }

267 268
private:

269

270 271 272
    ///////////////////////////////////////////////////////////////////////////////
    //                                  Computation
    ///////////////////////////////////////////////////////////////////////////////
273

274 275

    /** P2M computation
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
    * expansion_P2M_add
    * Multipole expansion with m charges q_i in Q_i=(rho_i, alpha_i, beta_i)
    *whose relative coordinates according to *p_center are:
    *Q_i - *p_center = (rho'_i, alpha'_i, beta'_i);
    *
    *For j=0..P, k=-j..j, we have:
    *
    *M_j^k = (-1)^j { sum{i=1..m} q_i Inner_j^k(rho'_i, alpha'_i, beta'_i) }
    *
    *However the extern loop is over the bodies (i=1..m) in our code and as an
    *intern loop we have: j=0..P, k=-j..j
    *
    *and the potential is then given by:
    *
    * Phi(x) = sum_{n=0}^{+} sum_{m=-n}^{n} M_n^m O_n^{-m} (x - *p_center)
    *
    */
293 294
    void particleToMultiPole(FComplex<FReal>* const cellMultiPole, const FPoint<FReal>& inPolePosition ,
                             const FPoint<FReal>& particlePosition, const FReal particlePhysicalValue){
berenger-bramas's avatar
berenger-bramas committed
295

296
        // Inner of Qi - Z0 => harmonic.result
297
        harmonic.computeInner( FSpherical<FReal>(particlePosition - inPolePosition) );
298

299
        FReal minus_one_pow_j = 1.0;    // (-1)^j => be in turn 1 and -1
300
        const FReal qParticle = particlePhysicalValue; // q in the formula
301
        int index_j_k = 0; // p_exp_term & p_Y_term
302

303 304 305 306 307 308
        // J from 0 to P
        for(int j = 0 ; j <= devP ; ++j){
            // k from 0 to J
            for(int k = 0 ; k <= j ; ++k, ++index_j_k){
                harmonic.result(index_j_k).mulRealAndImag( qParticle * minus_one_pow_j );
                cellMultiPole[index_j_k] += harmonic.result(index_j_k);
309
            }
310
            // (-1)^J => -1 becomes 1 or 1 becomes -1
311 312 313 314
            minus_one_pow_j = -minus_one_pow_j;
        }
    }

315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
    /* M2M
    *We compute the translation of multipole_exp_src from *p_center_of_exp_src to
    *p_center_of_exp_target, and add the result to multipole_exp_target.
    *
    * O_n^l (with n=0..P, l=-n..n) being the former multipole expansion terms
    * (whose center is *p_center_of_multipole_exp_src) we have for the new multipole
    * expansion terms (whose center is *p_center_of_multipole_exp_target):

    * M_j^k = sum{n=0..j}
    * sum{l=-n..n, |k-l|<=j-n}
    * O_n^l Inner_{j-n}^{k-l}(rho, alpha, beta)
    *
    * where (rho, alpha, beta) are the spherical coordinates of the vector :
    * p_center_of_multipole_exp_target - *p_center_of_multipole_exp_src
    *
    * Warning: if j-n < |k-l| we do nothing.
     */
332 333 334
    void multipoleToMultipole(FComplex<FReal>* const FRestrict multipole_exp_target,
                              const FComplex<FReal>* const FRestrict multipole_exp_src,
                              const FComplex<FReal>* const FRestrict M2M_Inner_transfer){
335

336 337
        // n from 0 to P
        for(int n = 0 ; n <= devP ; ++n ){
338
            // l<0 // (-1)^l
339
            FReal pow_of_minus_1_for_l = ( n & 1 ? FReal(-1.0) : FReal(1.0) );
340 341

            // O_n^l : here points on the source multipole expansion term of degree n and order |l|
342
            const int index_n = harmonic.getPreExpRedirJ(n);
343

344 345
            // l from -n to <0
            for(int l = -n ; l < 0 ; ++l){
346
                const FComplex<FReal> M_n__n_l = multipole_exp_src[index_n -l];
347

348 349
                // j from n to P
                for(int j = n ; j <= devP ; ++j ){
350
                    // M_j^k
351
                    const int index_j = harmonic.getPreExpRedirJ(j);
352
                    // Inner_{j-n}^{k-l} : here points on the M2M transfer function/expansion term of degree n-j and order |k-l|
353
                    const int index_j_n = harmonic.getPreExpRedirJ(j-n); /* k==0 */
354 355

                    // since n-j+l<0
356
                    for(int k = 0 ; k <= (j-n+l) ; ++k ){ // l<0 && k>=0 => k-l>0
357
                        const FComplex<FReal> I_j_n__k_l = M2M_Inner_transfer[index_j_n + k - l];
358 359

                        multipole_exp_target[index_j + k].incReal( pow_of_minus_1_for_l *
360 361
                                                    ((M_n__n_l.real() * I_j_n__k_l.real()) +
                                                     (M_n__n_l.imag() * I_j_n__k_l.imag())));
362
                        multipole_exp_target[index_j + k].incImag( pow_of_minus_1_for_l *
363 364
                                                    ((M_n__n_l.real() * I_j_n__k_l.imag()) -
                                                     (M_n__n_l.imag() * I_j_n__k_l.real())));
365 366 367

                     } // for k
                } // for j
368 369

                pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
370 371
            } // for l

372 373
            // l from 0 to n
            for(int l = 0 ; l <= n ; ++l){
374
                const FComplex<FReal> M_n__n_l = multipole_exp_src[index_n + l];
375

376 377 378
                // j from n to P
                for( int j = n ; j <= devP ; ++j ){
                    const int first_k = FMath::Max(0,n-j+l);
379
                    // (-1)^k
380
                    FReal pow_of_minus_1_for_k = static_cast<FReal>( first_k&1 ? -1.0 : 1.0 );
381
                    // M_j^k
382
                    const int index_j = harmonic.getPreExpRedirJ(j);
383
                    // Inner_{j-n}^{k-l} : here points on the M2M transfer function/expansion term of degree n-j and order |k-l|
384 385 386 387
                    const int index_j_n = harmonic.getPreExpRedirJ(j-n);

                    int k = first_k;
                    for(; k <= (j-n+l) && k < l ; ++k){ /* l>=0 && k-l<0 */
388
                        const FComplex<FReal> I_j_n__l_k = M2M_Inner_transfer[index_j_n + l - k];
389 390

                        multipole_exp_target[index_j + k].incReal( pow_of_minus_1_for_k * pow_of_minus_1_for_l *
391 392
                                                    ((M_n__n_l.real() * I_j_n__l_k.real()) +
                                                     (M_n__n_l.imag() * I_j_n__l_k.imag())));
393
                        multipole_exp_target[index_j + k].incImag(pow_of_minus_1_for_k * pow_of_minus_1_for_l *
394 395
                                                   ((M_n__n_l.imag() * I_j_n__l_k.real()) -
                                                    (M_n__n_l.real() * I_j_n__l_k.imag())));
396 397

                        pow_of_minus_1_for_k = -pow_of_minus_1_for_k;
398 399
                    } // for k

400
                    for(/* k = l */; k <= (j - n + l) ; ++k){ // l>=0 && k-l>=0
401
                        const FComplex<FReal> I_j_n__k_l = M2M_Inner_transfer[index_j_n + k - l];
402 403

                        multipole_exp_target[index_j + k].incReal(
404 405
                                (M_n__n_l.real() * I_j_n__k_l.real()) -
                                (M_n__n_l.imag() * I_j_n__k_l.imag()));
406
                        multipole_exp_target[index_j + k].incImag(
407 408
                                (M_n__n_l.imag() * I_j_n__k_l.real()) +
                                (M_n__n_l.real() * I_j_n__k_l.imag()));
berenger-bramas's avatar
berenger-bramas committed
409

410 411
                    } // for k
                } // for j
412 413

                pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
414 415 416 417 418 419
            } // for l
        } // for n
    }


    /** L2L
420
      *We compute the shift of local_exp_src from *p_center_of_exp_src to
421
      *p_center_of_exp_target, and set the result to local_expansion_target.
422 423 424 425 426 427 428 429 430 431 432 433 434
      *
      *O_n^l (with n=0..P, l=-n..n) being the former local expansion terms
      *(whose center is *p_center_of_exp_src) we have for the new local
      *expansion terms (whose center is *p_center_of_exp_target):
      *
      *L_j^k = sum{n=j..P}
      *sum{l=-n..n}
      *O_n^l Inner_{n-j}^{l-k}(rho, alpha, beta)
      *
      *where (rho, alpha, beta) are the spherical coordinates of the vector :
      *p_center_of_exp_target - *p_center_of_exp_src
      *
      *Warning: if |l-k| > n-j, we do nothing.
435
      */
436
    void localToLocal(FComplex<FReal>* const FRestrict local_expansion_target, const FComplex<FReal>* const FRestrict local_exp_src,
437
                      const FComplex<FReal>* const FRestrict L2L_tranfer){
438
        // L_j^k
439 440 441
        int index_j_k = 0;

        for (int j = 0 ; j <= devP ; ++j ){
442 443
            // (-1)^k
            FReal pow_of_minus_1_for_k = 1.0;
444 445

            for (int k = 0 ; k <= j ; ++k, ++index_j_k ){
446
                FComplex<FReal> L_j_k = local_expansion_target[index_j_k];
447 448

                for (int n=j; n <= devP;++n){
449
                    // O_n^l : here points on the source multipole expansion term of degree n and order |l|
450
                    const int index_n = harmonic.getPreExpRedirJ(n);
451

452
                    int l = n - j + k;
453
                    // Inner_{n-j}^{l-k} : here points on the L2L transfer function/expansion term of degree n-j and order |l-k|
454
                    const int index_n_j = harmonic.getPreExpRedirJ(n-j);
455

456
                    for(/*l = n - j + k*/ ; l-k > 0 ;  --l){ /* l>0 && l-k>0 */
457 458
                        const FComplex<FReal> L_j_l = local_exp_src[index_n + l];
                        const FComplex<FReal> I_l_j__l_k = L2L_tranfer[index_n_j  + l - k];
459

460 461 462 463
                        L_j_k.incReal( (L_j_l.real() * I_l_j__l_k.real()) -
                                                    (L_j_l.imag() * I_l_j__l_k.imag()));
                        L_j_k.incImag( (L_j_l.imag() * I_l_j__l_k.real()) +
                                                    (L_j_l.real() * I_l_j__l_k.imag()));
berenger-bramas's avatar
berenger-bramas committed
464

465 466 467
                    }

                    // (-1)^l
468 469
                    FReal pow_of_minus_1_for_l = ((l&1) ? FReal(-1.0) : FReal(1.0));
                    for(/*l = k*/; l>0 && l>=j-n+k; --l){ /* l>0 && l-k<=0 */
470 471
                        const FComplex<FReal> L_j_l = local_exp_src[index_n + l];
                        const FComplex<FReal> I_l_j__l_k = L2L_tranfer[index_n_j  - l + k];
472 473

                        L_j_k.incReal( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
474 475
                                                    ((L_j_l.real() * I_l_j__l_k.real()) +
                                                     (L_j_l.imag() * I_l_j__l_k.imag())));
476
                        L_j_k.incImag( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
477 478
                                                    ((L_j_l.imag() * I_l_j__l_k.real()) -
                                                     (L_j_l.real() * I_l_j__l_k.imag())));
479 480

                        pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
481 482 483
                     }

                    // l<=0 && l-k<=0
484
                    for(/*l = 0 ou l = j-n+k-1*/; l>=j-n+k; --l){
485 486
                        const FComplex<FReal> L_j_l = local_exp_src[index_n - l];
                        const FComplex<FReal> I_l_j__l_k = L2L_tranfer[index_n_j  - l + k];
487 488

                        L_j_k.incReal( pow_of_minus_1_for_k *
489 490
                                                    ((L_j_l.real() * I_l_j__l_k.real()) -
                                                     (L_j_l.imag() * I_l_j__l_k.imag())));
491
                        L_j_k.decImag( pow_of_minus_1_for_k *
492 493
                                                    ((L_j_l.imag() * I_l_j__l_k.real()) +
                                                     (L_j_l.real() * I_l_j__l_k.imag())));
berenger-bramas's avatar
berenger-bramas committed
494

495

496
                    }
497 498
                }//n

499
                local_expansion_target[index_j_k] = L_j_k;
500 501 502 503

                pow_of_minus_1_for_k = -pow_of_minus_1_for_k;
            }//k
        }//j
504 505 506 507
    }

    /** L2P
      */
508 509
    void localToParticle(const FPoint<FReal>& local_position,const FComplex<FReal>*const local_exp,
                         const FPoint<FReal>& particlePosition,
510 511
                         const FReal physicalValue, FReal*const potential,
                         FReal*const forcesX,FReal*const forcesY,FReal*const forcesZ){
512 513
        //--------------- Forces ----------------//

514 515 516
        FReal force_vector_in_local_base_x = 0.0;
        FReal force_vector_in_local_base_y = 0.0;
        FReal force_vector_in_local_base_z = 0.0;
517

518
        const FSpherical<FReal> spherical(particlePosition - local_position);
519
        harmonic.computeInnerTheta( spherical );
520 521 522 523 524
//        std::cout << "  L2P:"<<std::endl
//        		  << "        Centre: " << 		local_position <<std::endl
//        		  << "        PArt: " << 		particlePosition <<std::endl
//        		  << "        Diff           " << 		particlePosition - local_position<<std::endl
//        		  << "        Spherical Diff " << 		spherical <<std::endl ;
525

berenger-bramas's avatar
berenger-bramas committed
526
        int index_j_k = 1;
527

berenger-bramas's avatar
berenger-bramas committed
528
        for (int j = 1 ; j <= devP ; ++j ){
529 530
            {
                // k=0:
531
                // F_r:
532
                const FReal exp_term_aux_real = ( (harmonic.result(index_j_k).real() * local_exp[index_j_k].real()) - (harmonic.result(index_j_k).imag() * local_exp[index_j_k].imag()) );
533 534 535 536 537
                //const FReal exp_term_aux_imag = ( (harmonic.result(index_j_k).getReal() * local_exp[index_j_k].getImag()) + harmonic.result(index_j_k).getImag() * local_exp[index_j_k].getReal()) );
                force_vector_in_local_base_x = ( force_vector_in_local_base_x  + FReal(j) * exp_term_aux_real );
            }
            {
                // F_phi: k=0 => nothing to do for F_phi
538
                // F_theta:
539
                const FReal exp_term_aux_real = ( (harmonic.resultThetaDerivated(index_j_k).real() * local_exp[index_j_k].real()) - (harmonic.resultThetaDerivated(index_j_k).imag() * local_exp[index_j_k].imag()) );
540 541 542
                //const FReal exp_term_aux_imag = ( (harmonic.resultThetaDerivated(index_j_k).getReal() * local_exp[index_j_k].getImag()) + (harmonic.resultThetaDerivated(index_j_k).getImag() * local_exp[index_j_k].getReal()) );
                force_vector_in_local_base_y = ( force_vector_in_local_base_y + exp_term_aux_real );
            }
berenger-bramas's avatar
berenger-bramas committed
543

544
            ++index_j_k;
545

546 547 548 549
            // k>0:
            for (int k=1; k<=j ;++k, ++index_j_k){
                {
                    // F_r:
550 551
                    const FReal exp_term_aux_real = ( (harmonic.result(index_j_k).real() * local_exp[index_j_k].real()) - (harmonic.result(index_j_k).imag() * local_exp[index_j_k].imag()) );
                    const FReal exp_term_aux_imag = ( (harmonic.result(index_j_k).real() * local_exp[index_j_k].imag()) + (harmonic.result(index_j_k).imag() * local_exp[index_j_k].real()) );
552 553 554 555 556 557
                    force_vector_in_local_base_x = (force_vector_in_local_base_x  + FReal(2 * j) * exp_term_aux_real );
                    // F_phi:
                    force_vector_in_local_base_z = ( force_vector_in_local_base_z - FReal(2 * k) * exp_term_aux_imag);
                }
                {
                    // F_theta:
558
                    const FReal exp_term_aux_real = ( (harmonic.resultThetaDerivated(index_j_k).real() * local_exp[index_j_k].real()) - (harmonic.resultThetaDerivated(index_j_k).imag() * local_exp[index_j_k].imag()) );
559 560 561
                    //const FReal exp_term_aux_imag = ( (harmonic.resultThetaDerivated(index_j_k).getReal() * local_exp[index_j_k].getImag()) + (harmonic.resultThetaDerivated(index_j_k).getImag() * local_exp[index_j_k].getReal()) );
                    force_vector_in_local_base_y = (force_vector_in_local_base_y + FReal(2.0) * exp_term_aux_real );
                }
berenger-bramas's avatar
berenger-bramas committed
562

563
            }
564 565 566 567

        }
        // We want: - gradient(POTENTIAL_SIGN potential).
        // The -(- 1.0) computing is not the most efficient programming ...
berenger-bramas's avatar
berenger-bramas committed
568
        const FReal signe = 1.0;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
569
        //if( FMath::Epsilon < spherical.getR()){
570 571 572
        // The derivative wrt r is equivalent to use classical outer expansion multiplied by l and the result is divided by r
            force_vector_in_local_base_x = ( force_vector_in_local_base_x  * signe / spherical.getR());
            // classical definition of the derivatives
573 574
            force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR());
            force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta()));
BRAMAS Berenger's avatar
BRAMAS Berenger committed
575
        //}
576
        /////////////////////////////////////////////////////////////////////
577 578
//
            //spherical_position_Set_ph
579
        FReal ph = 	spherical.getPhiZero2Pi();
580 581 582
//        FReal ph = FMath::Fmod(spherical.getPhi(), FReal(2)*FMath::FPi<FReal>());
//        if (ph > M_PI) ph -= FReal(2) * FMath::FPi<FReal>();
//        if (ph < -M_PI + FMath::Epsilon)  ph += FReal(2) * FMath::FPi<FReal>();
583 584
//
//        //spherical_position_Set_th
585
//        FReal th = FMath::Fmod(spherical.getTheta(), FReal(2) * FMath::FPi<FReal>());
586
//        //FReal th = spherical.getTheta();
587 588 589
//        if (th < 0.0) th += 2*FMath::FPi<FReal>();
//        if (th > FMath::FPi<FReal>()){
//            th = 2*FMath::FPi<FReal>() - th;
590
//            //spherical_position_Set_ph(p, spherical_position_Get_ph(p) + M_PI);
591 592 593
//            ph = FMath::Fmod(ph + FMath::FPi<FReal>(), 2*FMath::FPi<FReal>());
//            if (ph > M_PI) ph -= 2*FMath::FPi<FReal>();
//            if (ph < -M_PI + FMath::Epsilon)  ph += 2 * FMath::FPi<FReal>();
594 595 596
//        }
//        //spherical_position_Set_r
//        //FReal rh = spherical.r;
BRAMAS Berenger's avatar
BRAMAS Berenger committed
597
        FAssertLF(spherical.getR() >= 0 , "R should be < 0!");
598

599 600
        const FReal cos_theta   = spherical.getCosTheta(); //FMath::Cos(th);
    //    const FReal cos_theta   =  FMath::Cos(th);
601
        const FReal cos_phi     = FMath::Cos(ph);
602 603 604
        const FReal sin_theta   = spherical.getSinTheta() ;// FMath::Sin(th);
   //     const FReal sin_theta   =  FMath::Sin(th);
       const FReal sin_phi     = FMath::Sin(ph);
605 606 607 608
//
        //
        // Formulae below are OK
        //
609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626
        FReal force_vector_tmp_x = (
                cos_phi * sin_theta * force_vector_in_local_base_x  +
                cos_phi * cos_theta * force_vector_in_local_base_y +
                (-sin_phi) * force_vector_in_local_base_z);

        FReal force_vector_tmp_y = (
                sin_phi * sin_theta * force_vector_in_local_base_x  +
                sin_phi * cos_theta * force_vector_in_local_base_y +
                cos_phi * force_vector_in_local_base_z);

        FReal force_vector_tmp_z = (
                cos_theta * force_vector_in_local_base_x +
                (-sin_theta) * force_vector_in_local_base_y);

        force_vector_tmp_x *= physicalValue;
        force_vector_tmp_y *= physicalValue;
        force_vector_tmp_z *= physicalValue;

627 628 629
        (*forcesX) += force_vector_tmp_x;
        (*forcesY) += force_vector_tmp_y;
        (*forcesZ) += force_vector_tmp_z;
630
        //
631
        //--------------- Potential ----------------//
632
        //
633
        FReal result = 0.0;
634
        index_j_k    = 0;
635 636

        for(int j = 0 ; j<= devP ; ++j){
637
            // k=0
638
            harmonic.result(index_j_k) *= local_exp[index_j_k];
639
            result += harmonic.result(index_j_k).real();
berenger-bramas's avatar
berenger-bramas committed
640

641
            ++index_j_k;
642 643

            // k>0
644 645
            for (int k = 1 ; k <= j ; ++k, ++index_j_k){
                harmonic.result(index_j_k) *= local_exp[index_j_k];
646
                result += 2 * harmonic.result(index_j_k).real();
647 648 649
            }
        }

650
        (*potential) += (result /* * physicalValue*/);
651 652
    }

653
public:
654 655 656
    /** Update a velocity of a particle
      *
      */
657 658 659 660 661
    void computeVelocity(ContainerClass*const FRestrict inParticles, const FReal DT){
        const FReal*const physicalValues = inParticles->getPhysicalValues();
        FReal*const forcesX = inParticles->getForcesX();
        FReal*const forcesY = inParticles->getForcesY();
        FReal*const forcesZ = inParticles->getForcesZ();
662
        FVector<FPoint<FReal>>& velocities = inParticles->getVelocities();
663

664
        for(FSize idxPart = 0 ; idxPart < inParticles->getNbParticles() ; ++idxPart){
665 666 667 668 669
            const FReal physicalValue = physicalValues[idxPart];
            // Coef = 1/m * time/2
            const FReal coef = (FReal(1.0)/physicalValue) * (DT/FReal(2.0));

            // velocity = velocity + forces * coef
670
            FPoint<FReal> forces_coef(forcesX[idxPart], forcesY[idxPart], forcesZ[idxPart]);
671 672 673
            forces_coef *= coef;
            velocities[idxPart] += (forces_coef);
        }
674 675 676 677 678
    }

    /** Update a position of a particle
      *
      */
679
    void updatePosition(ContainerClass*const FRestrict inParticles, const FReal DT){
680 681 682
        FReal*const positionsX = inParticles->getPositions()[0];
        FReal*const positionsY = inParticles->getPositions()[1];
        FReal*const positionsZ = inParticles->getPositions()[2];
683
        FVector<FPoint<FReal>>& velocities = inParticles->getVelocities();
684

685
        for(FSize idxPart = 0 ; idxPart < inParticles->getNbParticles() ; ++idxPart){
686
            FPoint<FReal> velocity_dt( velocities[idxPart] );
687 688 689 690 691
            velocity_dt *= DT;
            positionsX[idxPart] += velocity_dt.getX();
            positionsY[idxPart] += velocity_dt.getY();
            positionsZ[idxPart] += velocity_dt.getZ();
        }
692
    }
693 694 695
};


696
#endif //FABSTRACTSPHERICALKERNEL_HPP