FAdaptiveKernelWrapper.hpp 17.4 KB
Newer Older
1 2
#ifndef FAdaptiveKERNELWRAPPER_HPP
#define FAdaptiveKERNELWRAPPER_HPP
3

4 5
#include "Components/FAbstractKernels.hpp"
#include "Containers/FVector.hpp"
6

7
#include "FAdaptiveCell.hpp"
8 9 10 11



/**
12
 * This class is a wrapper to use the usual algorithm but in a Adaptive way.
13 14 15 16
 * The real computational kernel should be given in template.
 * It must propose the usual FMM Kernel operators and also the FAbstractKernels operators.
 *
 */
17 18
template <class RealAdaptiveKernel, class CellClass, class ContainerClass >
class FAdaptiveKernelWrapper : public FAbstractKernels<FAdaptiveCell<CellClass, ContainerClass>, ContainerClass>{
19
protected:
20
    RealAdaptiveKernel kernel;
21 22
public:
    template<typename... KernelParams>
23
    FAdaptiveKernelWrapper(KernelParams... parameters) : kernel(parameters...){
24 25
    }

26
    RealAdaptiveKernel& getKernel(){
27 28 29
        return kernel;
    }

30
    const RealAdaptiveKernel& getKernel() const{
31 32 33 34 35 36 37
        return kernel;
    }

    /** P2M is performed only if the kernel says that is it better than waiting for future development
      * A leaf cell contains the particles container if the development is not made,
      * else it has a multipole component.
      */
38 39 40
    void P2M(FAdaptiveCell<CellClass, ContainerClass>* const pole, const ContainerClass* const particles)  override {
        // Leaf is always Adaptive
        pole->setAdaptive(true);
41 42 43 44 45 46 47 48 49 50 51 52 53
        if( kernel.preferP2M(particles) ){
            // If it is better to compute the P2M at this level
            pole->setHaveDevelopment(true);
            kernel.P2M(pole->getRealCell(), particles);
        }
        else{
            // Else simply keep the current leaf
            pole->setHaveDevelopment(false);
            pole->addSubLeaf(particles);
        }
    }

    /** The M2M need to manage all the possibilities and to propagate the data from the leaf.
54
      * If there is one child, we are not Adaptive and keep the lower adaptive cell as target
55 56 57 58 59 60
      * If at least one child has a development (at any level) the current cell store the development of all the children
      * using M2M or P2M from the leaves.
      * If no children have developments we need to test if it is better to merge them in a development
      * or to store all the particles in order to perform the P2M later.
      * Finally if all the children have developments at level just under we can perform normal M2M
      */
61 62
    void M2M(FAdaptiveCell<CellClass, ContainerClass>* const FRestrict pole,
             const FAdaptiveCell<CellClass, ContainerClass>*const FRestrict *const FRestrict child, const int inLevel)  override {
63 64 65 66 67 68 69 70 71 72
        int nbChild       = 0;
        int lastChild     = 0;
        bool onlyParticlesCells         = true;
        FVector<const ContainerClass*> subLeaves;

        // Test all the children
        for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
            if(child[idxChild]){
                nbChild  += 1;
                lastChild = idxChild;
73 74
                // We aggregate all the leaves from the current child or its Adaptive cell
                if(child[idxChild]->isAdaptive() && !child[idxChild]->hasDevelopment()){
75 76
                    subLeaves.memocopy( child[idxChild]->getSubLeaves(), child[idxChild]->getNbSubLeaves());
                }
77 78 79
                else if(!child[idxChild]->isAdaptive() && !child[idxChild]->getSubAdaptiveCell()->hasDevelopment()){
                    subLeaves.memocopy( child[idxChild]->getSubAdaptiveCell()->getSubLeaves(),
                                        child[idxChild]->getSubAdaptiveCell()->getNbSubLeaves());
80 81 82 83 84 85 86
                }
                else{
                    // If a child is made of development
                    onlyParticlesCells = false;
                }
            }
        }
87
        // We need to aggregate if there are only particles and if the kernel says so
88 89
        const bool continueToAgregate = (onlyParticlesCells && kernel.preferP2M(inLevel, subLeaves.data(), subLeaves.getSize()));
        if(nbChild == 1){
90 91
            // One child means that the cell is not Adaptive
            pole->setAdaptive(false);
92
            pole->setHaveDevelopment(false);
93 94
            if(child[lastChild]->isAdaptive()){
                pole->setSubAdaptiveCell(child[lastChild], inLevel + 1);
95 96
            }
            else{
97 98
                pole->setSubAdaptiveCell(child[lastChild]->getSubAdaptiveCell(),
                                           child[lastChild]->getSubAdaptiveLevel());
99 100 101 102
            }
        }
        else if(onlyParticlesCells && continueToAgregate){
            // There are only particles and no development to do
103
            pole->setAdaptive(true);
104 105 106 107 108 109 110 111 112
            pole->setHaveDevelopment(false);
            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
                if(child[idxChild]){
                    pole->addSubLeaves(child[idxChild]->getSubLeaves(), child[idxChild]->getNbSubLeaves());
                }
            }
        }
        else{
            // There development to do from developments or particles
113
            pole->setAdaptive(true);
114 115 116 117 118 119 120
            pole->setHaveDevelopment(true);

            const CellClass* realChild[8] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
            int counterRealLowerCell      = 0;
            // Test each child
            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
                if(child[idxChild]){
121
                    if(child[idxChild]->isAdaptive()){
122
                        if(child[idxChild]->hasDevelopment()){
123
                            // If it is Adaptive and has development than we compute is using usual M2M
124 125 126 127
                            realChild[idxChild] = child[idxChild]->getRealCell();
                            counterRealLowerCell += 1;
                        }
                        else{
128
                            // If it is Adaptive and has not development than we compute is using P2M
129 130 131 132 133 134
                            for(int idxLeaf = 0 ; idxLeaf < child[idxChild]->getNbSubLeaves() ; ++idxLeaf){
                                kernel.P2M(pole->getRealCell(), inLevel, child[idxChild]->getSubLeaf(idxLeaf));
                            }
                        }
                    }
                    else{
135 136 137 138
                        // Else we need the Adaptive cell
                        const FAdaptiveCell<CellClass, ContainerClass>* lowerAdaptiveCell = child[idxChild]->getSubAdaptiveCell();
                        const int lowerAdaptiveLevel = child[idxChild]->getSubAdaptiveLevel();
                        if(lowerAdaptiveCell->hasDevelopment()){
139
                            // If it has development we perform a M2M
140 141
                            kernel.M2M(pole->getRealCell(), inLevel, lowerAdaptiveCell->getRealCell(),
                                       lowerAdaptiveLevel);
142 143 144 145
                        }
                        else{
                            // Else we perform P2M
                            for(int idxLeaf = 0 ; idxLeaf < child[idxChild]->getNbSubLeaves() ; ++idxLeaf){
146
                                kernel.P2M(pole->getRealCell(), inLevel, lowerAdaptiveCell->getSubLeaf(idxLeaf));
147 148 149 150 151 152 153 154 155 156 157 158 159
                            }
                        }
                    }
                }
            }
            // If there are usual M2M to do
            if( counterRealLowerCell ){
                kernel.M2M(pole->getRealCell(), realChild, inLevel);
            }
        }
    }


160 161
    /** The M2L should take into account if the current cell is Adaptive or not.
      * Else we should work on the sub Adaptive.
162 163 164
      * If it is composed or particles we have to perform M2P or P2P
      * Else we have to perform M2L or P2L.
      */
165 166
    void M2L(FAdaptiveCell<CellClass, ContainerClass>* const FRestrict local,
                const FAdaptiveCell<CellClass, ContainerClass>* distantNeighbors[343],
167 168 169 170
                const int /*size*/, const int inLevel)  override {
        // In case usual M2L can be done
        const CellClass* normalDistantNeighbors[343];
        int normalSize = 0;
171 172 173 174 175 176 177
        // The current Adaptive cell
        FAdaptiveCell<CellClass, ContainerClass>* currentAdaptiveCell = nullptr;
        int currentAdaptiveLevel       = -1;
        // If the current Adaptive cell is the current cell
        if(local->isAdaptive()){
            currentAdaptiveCell  = local;
            currentAdaptiveLevel = inLevel;
178 179 180 181 182
            // Then we may have some M2L to do
            memset(normalDistantNeighbors, 0, 343*sizeof(CellClass*));
        }
        else{
            // Else we are working with a lower cell
183 184
            currentAdaptiveCell = local->getSubAdaptiveCell();
            currentAdaptiveLevel= local->getSubAdaptiveLevel();
185 186 187 188
        }

        for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
            if(distantNeighbors[idxNeigh]){
189 190 191
                // If the current cell is Adaptive and the neighbor too
                if(distantNeighbors[idxNeigh]->isAdaptive() && local->isAdaptive()){
                    if(distantNeighbors[idxNeigh]->hasDevelopment() && currentAdaptiveCell->hasDevelopment()){
192 193 194 195
                        // If both have development than we can use usual M2L
                        normalDistantNeighbors[idxNeigh] = distantNeighbors[idxNeigh]->getRealCell();
                        normalSize += 1;
                    }
196
                    else if(currentAdaptiveCell->hasDevelopment()){
197 198
                        // If only current cell has development the neighbor has particles
                        for(int idxLeafSrc = 0 ; idxLeafSrc < distantNeighbors[idxNeigh]->getNbSubLeaves() ; ++idxLeafSrc){
199
                            kernel.P2L(currentAdaptiveCell->getRealCell(), currentAdaptiveLevel, distantNeighbors[idxNeigh]->getSubLeaf(idxLeafSrc));
200 201 202 203
                        }
                    }
                    else if(distantNeighbors[idxNeigh]->hasDevelopment()){
                        // If only current cell has particles the neighbor has development
204 205
                        for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
                            kernel.M2P(distantNeighbors[idxNeigh]->getRealCell(), currentAdaptiveLevel, currentAdaptiveCell->getSubLeaf(idxLeafTgt));
206 207 208 209
                        }
                    }
                    else{
                        // If both have particles
210
                        for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
211
                            for(int idxLeafSrc = 0 ; idxLeafSrc < distantNeighbors[idxNeigh]->getNbSubLeaves() ; ++idxLeafSrc){
212
                                kernel.P2P(currentAdaptiveCell->getSubLeaf(idxLeafTgt), distantNeighbors[idxNeigh]->getSubLeaf(idxLeafSrc));
213 214 215 216 217
                            }
                        }
                    }
                }
                else{
218 219 220 221 222 223
                    const FAdaptiveCell<CellClass, ContainerClass>* lowerAdaptiveCell = distantNeighbors[idxNeigh];
                    int lowerAdaptiveLevel       = inLevel;
                    // If we need to look at lower level to find the Adaptive cell
                    if(!distantNeighbors[idxNeigh]->isAdaptive()){
                        lowerAdaptiveCell  = distantNeighbors[idxNeigh]->getSubAdaptiveCell();
                        lowerAdaptiveLevel = distantNeighbors[idxNeigh]->getSubAdaptiveLevel();
224 225
                    }

226
                    if(lowerAdaptiveCell->hasDevelopment() && currentAdaptiveCell->hasDevelopment()){
227
                        // We are doing a M2L with distant interaction
228 229
                        kernel.M2L(currentAdaptiveCell->getRealCell(), currentAdaptiveLevel,
                            lowerAdaptiveCell->getRealCell(), lowerAdaptiveLevel);
230
                    }
231
                    else if(currentAdaptiveCell->hasDevelopment()){
232
                        // If only current cell has development the neighbor has particles
233 234
                        for(int idxLeafSrc = 0 ; idxLeafSrc < lowerAdaptiveCell->getNbSubLeaves() ; ++idxLeafSrc){
                            kernel.P2L(currentAdaptiveCell->getRealCell(), currentAdaptiveLevel, lowerAdaptiveCell->getSubLeaf(idxLeafSrc));
235 236
                        }
                    }
237
                    else if(lowerAdaptiveCell->hasDevelopment()){
238
                        // If only current cell has particles the neighbor has development
239 240
                        for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
                            kernel.M2P(lowerAdaptiveCell->getRealCell(), currentAdaptiveLevel, currentAdaptiveCell->getSubLeaf(idxLeafTgt));
241 242 243 244
                        }
                    }
                    else{
                        // If both have particles
245 246 247
                        for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
                            for(int idxLeafSrc = 0 ; idxLeafSrc < lowerAdaptiveCell->getNbSubLeaves() ; ++idxLeafSrc){
                                kernel.P2P(currentAdaptiveCell->getSubLeaf(idxLeafTgt), lowerAdaptiveCell->getSubLeaf(idxLeafSrc));
248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263
                            }
                        }
                    }
                }
            }
        }

        // If we work on the current cell and it has development
        if(normalSize){
            kernel.M2L(local->getRealCell(), normalDistantNeighbors, normalSize, inLevel);
        }
    }


    /** Nothing special */
    void finishedLevelM2L(const int level){
264
    		kernel.finishedLevelM2L(level);
265 266 267
    }


268 269 270
    /** If the current cell is not Adaptive we have nothing to do.
      * If it is Adaptive we have to test the children.
      * If a child is Adaptive than
271 272
      *     it has development and it is a normal L2L that should be performed
      *     or it has a list of leaves (particles) and we need L2P
273
      * Else it points to the lower Adaptive cell and we need a L2L or L2P as in the previous case
274
      */
275 276
    void L2L(const FAdaptiveCell<CellClass, ContainerClass>* const FRestrict local,
             FAdaptiveCell<CellClass, ContainerClass>* FRestrict * const FRestrict child, const int inLevel)  override {
277
        // If there is something on this cell
278
        if(local->isAdaptive()){
279 280 281 282 283 284
            // We store the usual cell
            CellClass* realChild[8] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
            int counterRealChild    = 0;
            // For each child
            for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
                if(child[idxChild]){
285 286
                    // If it is Adaptive (it holds development or particles)
                    if(child[idxChild]->isAdaptive()){
287 288 289 290 291 292 293 294 295 296 297 298 299
                        if(child[idxChild]->hasDevelopment()){
                            // We need to perform a usual L2L on it
                            realChild[idxChild] = child[idxChild]->getRealCell();
                            counterRealChild   += 1;
                        }
                        else {
                            // We need to propagate on the particles
                            for(int idxLeafSrc = 0 ; idxLeafSrc < child[idxChild]->getNbSubLeaves() ; ++idxLeafSrc){
                                kernel.L2P(local->getRealCell(), inLevel, child[idxChild]->getSubLeaf(idxLeafSrc));
                            }
                        }
                    }
                    else{
300 301 302 303
                        // Get the lower Adaptive cell
                        FAdaptiveCell<CellClass, ContainerClass>* lowerAdaptiveCell = child[idxChild]->getSubAdaptiveCell();
                        const int lowerAdaptiveLevel = child[idxChild]->getSubAdaptiveLevel();
                        if(lowerAdaptiveCell->hasDevelopment()){
304
                            // If it has a development we do a L2L with more than 1 level difference
305
                            kernel.L2L(local->getRealCell(), inLevel, lowerAdaptiveCell->getRealCell(), lowerAdaptiveLevel);
306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
                        }
                        else{
                            // Else we propagate on the particles
                            for(int idxLeafSrc = 0 ; idxLeafSrc < child[idxChild]->getNbSubLeaves() ; ++idxLeafSrc){
                                kernel.L2P(local->getRealCell(), inLevel, child[idxChild]->getSubLeaf(idxLeafSrc));
                            }
                        }
                    }
                }
            }
            // Perform the usual L2L
            if(counterRealChild){
                kernel.L2L(local->getRealCell(), realChild, inLevel);
            }
        }
    }

    /** We do a Local to Particles only if the local (leaf) cell has some development */
324
    void L2P(const FAdaptiveCell<CellClass, ContainerClass>* const local, ContainerClass* const particles)  override {
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
        if(local->hasDevelopment()){
            kernel.L2P(local->getRealCell(), particles);
        }
    }

    /** This is a normal P2P */
    void P2P(const FTreeCoordinate& inLeafPosition,
             ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
             ContainerClass* const directNeighborsParticles[27], const int size)  override {
        kernel.P2P(inLeafPosition, targets, sources, directNeighborsParticles, size);
    }

    /** This is a normal P2P */
    void P2PRemote(const FTreeCoordinate& inLeafPosition,
                   ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                   ContainerClass* const directNeighborsParticles[27], const int size) override {
        kernel.P2PRemote(inLeafPosition, targets, sources, directNeighborsParticles, size);
    }

};

#endif