.. _program_listing_file_include_test.h:

Program Listing for File test.h
===============================

|exhale_lsh| :ref:`Return to documentation for file <file_include_test.h>` (``include/test.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #ifndef test_h
   #define test_h
   #include <cassert>
   #include "exafmm_t.h"
   #include "fmm_base.h"
   
   namespace exafmm_t {
     template <typename T>
     class DummyFmm : public FmmBase<T> {
     public:
       DummyFmm() {}
       DummyFmm(int ncrit_) { this->p = 1; this->nsurf = 1; this->ncrit = ncrit_;}
   
       ivec3 octant2coord(int octant) {
         ivec3 rel_coord;
         rel_coord[0] = octant & 1 ? 1 : -1;
         rel_coord[1] = octant & 2 ? 1 : -1;
         rel_coord[2] = octant & 4 ? 1 : -1;
         return rel_coord;
       } 
   
       void P2M(NodePtrs<T>& leafs) {
   #pragma omp parallel for
         for(size_t i=0; i<leafs.size(); ++i) {
           Node<T>* leaf = leafs[i];
           leaf->up_equiv[0] += leaf->nsrcs;
         }
       }
   
       void M2M(Node<T>* node) {
         if(node->is_leaf) return;
         for(int octant=0; octant<8; ++octant) {  
           if(node->children[octant])
   #pragma omp task untied
             M2M(node->children[octant]);
         }
   #pragma omp taskwait
         for(int octant=0; octant<8; octant++) {
           if(node->children[octant]) {
             Node<T>* child = node->children[octant];
             node->up_equiv[0] += child->up_equiv[0];
           }
         }
       }
   
       void M2L(NodePtrs<T>& nonleafs) {
   #pragma omp parallel for schedule(dynamic)
         for(size_t i=0; i<nonleafs.size(); ++i) {
           Node<T>* trg_parent = nonleafs[i];
           NodePtrs<T>& M2L_list = trg_parent->M2L_list;
           for(size_t j=0; j<M2L_list.size(); ++j) {
             Node<T>* src_parent = M2L_list[j];
             if(src_parent) {
               // find src_parent's relative position to trg_parent
               // NodePtrs<T>::iterator
               auto it = std::find(trg_parent->colleagues.begin(),
                                   trg_parent->colleagues.end(), src_parent);
               assert(it != trg_parent->colleagues.end());   // src_parent has to be trg_parent's colleague
               int colleague_idx = it - trg_parent->colleagues.begin();
               ivec3 parent_rel_coord;
               parent_rel_coord[0] = colleague_idx % 3 - 1;
               parent_rel_coord[1] = (colleague_idx/3) % 3 - 1;
               parent_rel_coord[2] = (colleague_idx/9) % 3 - 1;
               for(int src_octant=0; src_octant<8; ++src_octant) {
                 Node<T>* src_child = src_parent->children[src_octant];
                 ivec3 src_child_coord = octant2coord(src_octant);
                 for(int trg_octant=0; trg_octant<8; ++trg_octant) {
                   Node<T>* trg_child = trg_parent->children[trg_octant];
                   ivec3 trg_child_coord = octant2coord(trg_octant);
                   if(src_child && trg_child) {
                     ivec3 rel_coord = parent_rel_coord*2 + (src_child_coord - trg_child_coord) / 2;  // calculate relative coords between children
                     bool is_colleague = true;
                     for(int d=0; d<3; ++d) {
                       if(rel_coord[d]>1 || rel_coord[d]<-1)
                         is_colleague = false;
                     }
                     if(!is_colleague)  // perform M2L if they are not neighbors
                       trg_child->dn_equiv[0] += src_child->up_equiv[0];
                   }
                 }
               }
             }
           }
         }
       }
   
       void P2L(Nodes<T>& nodes) {
   #pragma omp parallel for schedule(dynamic)
         for(size_t i=0; i<nodes.size(); ++i) {
           NodePtrs<T>& P2L_list = nodes[i].P2L_list;
           for(size_t j=0; j<P2L_list.size(); ++j) {
             nodes[i].dn_equiv[0] += P2L_list[j]->nsrcs;
           }
         }
       }
   
       void M2P(NodePtrs<T>& leafs) {
   #pragma omp parallel for schedule(dynamic)
         for(size_t i=0; i<leafs.size(); ++i) {
           if (leafs[i]->ntrgs == 0) continue;
           NodePtrs<T>& M2P_list = leafs[i]->M2P_list;
           for (size_t j=0; j<M2P_list.size(); ++j) {
             leafs[i]->trg_value[0] += M2P_list[j]->up_equiv[0];
           }
         }
       }
   
       void L2L(Node<T>* node) {
         if(node->is_leaf) return;
         for(int octant=0; octant<8; octant++) {
           if(node->children[octant]) {
             Node<T>* child = node->children[octant];
             child->dn_equiv[0] += node->dn_equiv[0];
           }
         }
         for(int octant=0; octant<8; ++octant) {
           if(node->children[octant])
   #pragma omp task untied
             L2L(node->children[octant]);
         }
   #pragma omp taskwait
       }
   
       void L2P(NodePtrs<T>& leafs) {
   #pragma omp parallel for
         for(size_t i=0; i<leafs.size(); ++i) {
           Node<T>* leaf = leafs[i];
           if (leaf->ntrgs==0) continue;
           leaf->trg_value[0] += leaf->dn_equiv[0];
         }
       }
   
       void P2P(NodePtrs<T>& leafs) {
   #pragma omp parallel for schedule(dynamic)
         for(size_t i=0; i<leafs.size(); ++i) {
           Node<T>* leaf = leafs[i];
           if (leaf->ntrgs==0) continue;
           NodePtrs<T>& P2P_list = leaf->P2P_list;
           for(size_t j=0; j<P2P_list.size(); ++j) {
             leaf->trg_value[0] += P2P_list[j]->nsrcs;
           }
         }
       }
   
       // below are the virtual methods define in FmmBase class
       void potential_P2P(RealVec& src_coord, std::vector<T>& src_value,
                          RealVec& trg_coord, std::vector<T>& trg_value) {}
   
   
       void gradient_P2P(RealVec& src_coord, std::vector<T>& src_value,
                         RealVec& trg_coord, std::vector<T>& trg_value) {}
   
       void M2L(Nodes<T>& nodes) {}
     };
   
     template <typename T>
     void set_children(Node<T>* parent, Node<T>* first_child) {
       parent->is_leaf = false;
       for (int octant=0; octant<8; ++octant) {
         Node<T>* child = first_child + octant;
         child->octant = octant;
         child->parent = parent;
         child->level = parent->level + 1;
         child->x = parent->x;
         child->r = parent->r / 2;
         for (int d=0; d<3; d++) {
           child->x[d] += child->r * (((octant & 1 << d) >> d) * 2 - 1);
         }
         parent->children.push_back(child);
       }
     }
   }// end namespace
   #endif