FSphericalKernel.hpp 9.53 KB
Newer Older
1
// ===================================================================================
2 3 4 5 6 7 8 9
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
10 11 12 13
// ===================================================================================
#ifndef FSPHERICALKERNEL_HPP
#define FSPHERICALKERNEL_HPP

14
#include "FAbstractSphericalKernel.hpp"
15
#include "../../Utils/FMemUtils.hpp"
16 17 18

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
19
* This class is the basic spherical harmonic kernel
20 21
*/
template< class ParticleClass, class CellClass, class ContainerClass>
22 23 24
class FSphericalKernel : public FAbstractSphericalKernel<ParticleClass,CellClass,ContainerClass> {
protected:
    typedef FAbstractSphericalKernel<ParticleClass,CellClass,ContainerClass> Parent;
25

26
    const int devM2lP;               //< A secondary P
27

28
    FSmartPointer<FComplexe*> preM2LTransitions;   //< The pre-computation for the M2L based on the level and the 189 possibilities
29

30 31 32
    /** To access te pre computed M2L transfer vector */
    int indexM2LTransition(const int idxX,const int idxY,const int idxZ) const {
        return (((((idxX+3) * 7) + (idxY+3)) * 7 ) + (idxZ+3)) * devM2lP;
33 34 35 36 37 38
    }

    /** Alloc and init pre-vectors*/
    void allocAndInit(){
        // M2L transfer, there is a maximum of 3 neighbors in each direction,
        // so 6 in each dimension
39 40
        preM2LTransitions = new FComplexe*[Parent::treeHeight];
        memset(preM2LTransitions.getPtr(), 0, sizeof(FComplexe*) * (Parent::treeHeight));
41
        // We start from the higher level
42 43
        FReal treeWidthAtLevel = Parent::boxWidth;
        for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
44
            // Allocate data for this level
45
            preM2LTransitions[idxLevel] = new FComplexe[(7 * 7 * 7) * devM2lP];
46
            // Precompute transfer vector
47 48 49 50
            for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
                for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
                    for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
                        if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
51
                            const FPoint relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
52
                            Parent::harmonic.computeOuter(FSpherical(relativePos));
53
                            FMemUtils::copyall<FComplexe>(&preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)], Parent::harmonic.result(), Parent::harmonic.getExpSize());
54 55 56 57
                        }
                    }
                }
            }
58
            // We divide the bow per 2 when we go down
59 60 61 62 63 64
            treeWidthAtLevel /= 2;
        }
    }


public:
65 66 67 68 69 70
    /** Constructor
      * @param inDevP the polynomial degree
      * @param inThreeHeight the height of the tree
      * @param inBoxWidth the size of the simulation box
      * @param inPeriodicLevel the number of level upper to 0 that will be requiried
      */
71 72
    FSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter)
        : Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
73 74
          devM2lP(int(((inDevP*2)+1) * ((inDevP*2)+2) * 0.5)),
          preM2LTransitions(0) {
75 76 77 78 79
        allocAndInit();
    }

    /** Copy constructor */
    FSphericalKernel(const FSphericalKernel& other)
80 81 82
        : Parent(other), devM2lP(other.devM2lP),
          preM2LTransitions(other.preM2LTransitions) {

83 84
    }

85
    /** Destructor */
86
    ~FSphericalKernel(){
87
        if( preM2LTransitions.isLast() ){
88
            FMemUtils::DeleteAll(preM2LTransitions.getPtr(), Parent::treeHeight);
89
        }
90 91 92
    }

    /** M2L with a cell and all the existing neighbors */
93 94
    void M2L(CellClass* const FRestrict pole, const CellClass* distantNeighbors[343],
             const int /*size*/, const int inLevel) {
95
        // For all neighbors compute M2L
96 97
        for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
            if( distantNeighbors[idxNeigh] ){
98
                const FComplexe* const transitionVector = &preM2LTransitions[inLevel][idxNeigh * devM2lP];
99 100
                multipoleToLocal(pole->getLocal(), distantNeighbors[idxNeigh]->getMultipole(), transitionVector);
            }
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
        }
    }


    /** M2L
    *We compute the conversion of multipole_exp_src in *p_center_of_exp_src to
    *a local expansion in *p_center_of_exp_target, and add the result to local_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 local
    *expansion terms (whose center is *p_center_of_local_exp_target):
    *
    *L_j^k = sum{n=0..+}
    *sum{l=-n..n}
    *O_n^l Outer_{j+n}^{-k-l}(rho, alpha, beta)
    *
    *where (rho, alpha, beta) are the spherical coordinates of the vector :
    *p_center_of_local_exp_src - *p_center_of_multipole_exp_target
    *
    *Remark: here we have always j+n >= |-k-l|
      *
      */
    void multipoleToLocal(FComplexe*const FRestrict local_exp, const FComplexe* const FRestrict multipole_exp_src,
                          const FComplexe* const FRestrict M2L_Outer_transfer){
        int index_j_k = 0;

        // L_j^k
        // HPMSTART(51, "M2L computation (loops)");
        // j from 0 to P
130
        for (int j = 0 ; j <= Parent::devP ; ++j){
131 132 133 134 135 136 137 138 139 140
            // (-1)^k
            FReal pow_of_minus_1_for_k = 1.0;
            //k from 0 to j
            for (int k = 0 ; k <= j ; ++k, ++index_j_k){
                // (-1)^n
                FReal pow_of_minus_1_for_n = 1.0;

                // work with a local variable
                FComplexe L_j_k = local_exp[index_j_k];
                // n from 0 to P
141
                for (int n = 0 ; n <= /*devP or*/ Parent::devP-j ; ++n){
142
                    // O_n^l : here points on the source multipole expansion term of degree n and order |l|
143
                    const int index_n = Parent::harmonic.getPreExpRedirJ(n);
144 145

                    // Outer_{j+n}^{-k-l} : here points on the M2L transfer function/expansion term of degree j+n and order |-k-l|
146
                    const int index_n_j = Parent::harmonic.getPreExpRedirJ(n+j);
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205

                    FReal pow_of_minus_1_for_l = pow_of_minus_1_for_n; // (-1)^l

                    // We start with l=n (and not l=-n) so that we always set p_Outer_term to a correct value in the first loop.
                    int l = n;
                    for(/* l = n */ ; l > 0 ; --l){ // we have -k-l<0 and l>0
                        const FComplexe M_n_l = multipole_exp_src[index_n + l];
                        const FComplexe O_n_j__k_l = M2L_Outer_transfer[index_n_j + k + l];

                        L_j_k.incReal( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
                                                    ((M_n_l.getReal() * O_n_j__k_l.getReal()) +
                                                     (M_n_l.getImag() * O_n_j__k_l.getImag())));
                        L_j_k.incImag( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
                                                    ((M_n_l.getImag() * O_n_j__k_l.getReal()) -
                                                     (M_n_l.getReal() * O_n_j__k_l.getImag())));

                        pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
                    }

                    for(/* l = 0 */; l >= -n &&  (-k-l) < 0 ; --l){ // we have -k-l<0 and l<=0
                        const FComplexe M_n_l = multipole_exp_src[index_n - l];
                        const FComplexe O_n_j__k_l = M2L_Outer_transfer[index_n_j + k + l];

                        L_j_k.incReal( pow_of_minus_1_for_k *
                                                    ((M_n_l.getReal() * O_n_j__k_l.getReal()) -
                                                     (M_n_l.getImag() * O_n_j__k_l.getImag())));
                        L_j_k.decImag(  pow_of_minus_1_for_k *
                                                     ((M_n_l.getImag() * O_n_j__k_l.getReal()) +
                                                      (M_n_l.getReal() * O_n_j__k_l.getImag())));

                        pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
                    }

                    for(/*l = -n-1 or l = -k-1 */; l >= -n ; --l){ // we have -k-l>=0 and l<=0
                        const FComplexe M_n_l = multipole_exp_src[index_n - l];
                        const FComplexe O_n_j__k_l = M2L_Outer_transfer[index_n_j - (k + l)];

                        L_j_k.incReal( pow_of_minus_1_for_l *
                                                    ((M_n_l.getReal() * O_n_j__k_l.getReal()) +
                                                     (M_n_l.getImag() * O_n_j__k_l.getImag())));
                        L_j_k.incImag( pow_of_minus_1_for_l *
                                                    ((M_n_l.getReal() * O_n_j__k_l.getImag()) -
                                                     (M_n_l.getImag() * O_n_j__k_l.getReal())));

                        pow_of_minus_1_for_l = -pow_of_minus_1_for_l;
                    }

                    pow_of_minus_1_for_n = -pow_of_minus_1_for_n;
                }//n

                // put in the local vector
                local_exp[index_j_k] = L_j_k;

                pow_of_minus_1_for_k = -pow_of_minus_1_for_k;
            }//k
        }
    }
};

206
#endif // FSPHERICALKERNEL_HPP