FSphericalKernel.hpp 9.36 KB
Newer Older
1
// See LICENCE file at project root
2 3 4
#ifndef FSPHERICALKERNEL_HPP
#define FSPHERICALKERNEL_HPP

5
#include "FAbstractSphericalKernel.hpp"
6
#include "../../Utils/FMemUtils.hpp"
7 8 9

/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
10
* This class is the basic spherical harmonic kernel
11
*/
12 13
template< class FReal, class CellClass, class ContainerClass>
class FSphericalKernel : public FAbstractSphericalKernel<FReal, CellClass,ContainerClass> {
14
protected:
15
    typedef FAbstractSphericalKernel<FReal, CellClass,ContainerClass> Parent;
16

17
    const int devM2lP;               //< A secondary P
18

19
    FSmartPointer<FComplex<FReal>*> preM2LTransitions;   //< The pre-computation for the M2L based on the level and the 189 possibilities
20

21 22
    /** To access te pre computed M2L transfer vector */
    int indexM2LTransition(const int idxX,const int idxY,const int idxZ) const {
23
        return (((((idxX+3) * 7) + (idxY+3)) * 7 ) + (idxZ+3)) * devM2lP;
24 25 26 27
    }

    /** Alloc and init pre-vectors*/
    void allocAndInit(){
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
        // M2L transfer, there is a maximum of 3 neighbors in each direction,
        // so 6 in each dimension
        preM2LTransitions = new FComplex<FReal>*[Parent::treeHeight];
        memset(preM2LTransitions.getPtr(), 0, sizeof(FComplex<FReal>*) * (Parent::treeHeight));
        // We start from the higher level
        FReal treeWidthAtLevel = Parent::boxWidth;
        for(int idxLevel = 0 ; idxLevel < Parent::treeHeight ; ++idxLevel ){
            // Allocate data for this level
            preM2LTransitions[idxLevel] = new FComplex<FReal>[(7 * 7 * 7) * devM2lP];
            // Precompute transfer vector
            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){
                            const FPoint<FReal> relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
                            Parent::harmonic.computeOuter(FSpherical<FReal>(relativePos));
                            FMemUtils::copyall<FComplex<FReal>>(&preM2LTransitions[idxLevel][indexM2LTransition(idxX,idxY,idxZ)], Parent::harmonic.result(), Parent::harmonic.getExpSize());
                        }
                    }
                }
            }
            // We divide the bow per 2 when we go down
            treeWidthAtLevel /= 2;
        }
52 53 54 55
    }


public:
56 57 58 59 60 61
    /** 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
      */
62
    FSphericalKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter)
63 64 65 66
        : Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter),
          devM2lP(int(((inDevP*2)+1) * ((inDevP*2)+2) * 0.5)),
          preM2LTransitions(nullptr) {
        allocAndInit();
67 68 69 70
    }

    /** Copy constructor */
    FSphericalKernel(const FSphericalKernel& other)
71 72
        : Parent(other), devM2lP(other.devM2lP),
          preM2LTransitions(other.preM2LTransitions) {
73

74 75
    }

76
    /** Destructor */
77
    ~FSphericalKernel(){
78
        if( preM2LTransitions.isLast() ){
79
            FMemUtils::DeleteAllArray(preM2LTransitions.getPtr(), Parent::treeHeight);
80
        }
81 82 83
    }

    /** M2L with a cell and all the existing neighbors */
84 85 86 87 88 89 90
    template<class SymbolicData>
    void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
             const SymbolicData* const TargetSymb,
             const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
             const SymbolicData* const FRestrict /*SourceSymbs*/[],
             const int neighborPositions[],
             const int inSize)
91
    {
92
        // For all neighbors compute M2L
93
        int level = static_cast<int>(TargetSymb->getLevel());
94 95
        for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
            const int idxNeigh = neighborPositions[idxExistingNeigh];
96 97
            const FComplex<FReal>* const transitionVector = &preM2LTransitions[level][idxNeigh * devM2lP];
            multipoleToLocal(TargetExpansion->get(), SourceMultipoles[idxExistingNeigh]->get(), transitionVector);
98
        }
99 100 101 102 103
    }


    /** M2L
    *We compute the conversion of multipole_exp_src in *p_center_of_exp_src to
104
    *a local expansion in *p_center_of_exp_target, and add the result to local_expansion_target.
105 106 107
    *
    *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
108
    *expansion terms (whose center is *p_center_of_local_expansion_target):
109 110 111 112 113 114 115 116 117 118 119
    *
    *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|
      *
      */
120
    void multipoleToLocal(FComplex<FReal>*const FRestrict local_exp, const FComplex<FReal>* const FRestrict multipole_exp_src,
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
                          const FComplex<FReal>* const FRestrict M2L_Outer_transfer){
        int index_j_k = 0;

        // L_j^k
        // HPMSTART(51, "M2L computation (loops)");
        // j from 0 to P
        for (int j = 0 ; j <= Parent::devP ; ++j){
            // (-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
                FComplex<FReal> L_j_k = local_exp[index_j_k];
                // n from 0 to P or do P-j
                //for (int n = 0 ; n <= Parent::devP /*or*/ /*Parent::devP-j*/ ; ++n){
                for (int n = 0 ; n <= Parent::devP-j ; ++n){  // faster than double height Parent::devP
                    // O_n^l : here points on the source multipole expansion term of degree n and order |l|
                    const int index_n = Parent::harmonic.getPreExpRedirJ(n);

                    // Outer_{j+n}^{-k-l} : here points on the M2L transfer function/expansion term of degree j+n and order |-k-l|
                    const int index_n_j = Parent::harmonic.getPreExpRedirJ(n+j);

                    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 FComplex<FReal> M_n_l = multipole_exp_src[index_n + l];
                        const FComplex<FReal> 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 *
155 156
                                       ((M_n_l.real() * O_n_j__k_l.real()) +
                                        (M_n_l.imag() * O_n_j__k_l.imag())));
157
                        L_j_k.incImag( pow_of_minus_1_for_l * pow_of_minus_1_for_k *
158 159
                                       ((M_n_l.imag() * O_n_j__k_l.real()) -
                                        (M_n_l.real() * O_n_j__k_l.imag())));
160 161 162 163 164 165 166 167 168

                        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 FComplex<FReal> M_n_l = multipole_exp_src[index_n - l];
                        const FComplex<FReal> O_n_j__k_l = M2L_Outer_transfer[index_n_j + k + l];

                        L_j_k.incReal( pow_of_minus_1_for_k *
169 170
                                       ((M_n_l.real() * O_n_j__k_l.real()) -
                                        (M_n_l.imag() * O_n_j__k_l.imag())));
171
                        L_j_k.decImag(  pow_of_minus_1_for_k *
172 173
                                        ((M_n_l.imag() * O_n_j__k_l.real()) +
                                         (M_n_l.real() * O_n_j__k_l.imag())));
174 175 176 177 178 179 180 181 182

                        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 FComplex<FReal> M_n_l = multipole_exp_src[index_n - l];
                        const FComplex<FReal> O_n_j__k_l = M2L_Outer_transfer[index_n_j - (k + l)];

                        L_j_k.incReal( pow_of_minus_1_for_l *
183 184
                                       ((M_n_l.real() * O_n_j__k_l.real()) +
                                        (M_n_l.imag() * O_n_j__k_l.imag())));
185
                        L_j_k.incImag( pow_of_minus_1_for_l *
186 187
                                       ((M_n_l.real() * O_n_j__k_l.imag()) -
                                        (M_n_l.imag() * O_n_j__k_l.real())));
188 189 190 191 192 193 194 195 196 197 198 199 200

                        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
        }
201 202 203
    }
};

204
#endif // FSPHERICALKERNEL_HPP