GenericTreeNode.h

Go to the documentation of this file.

 
#ifndef GENERICTREENODE_H
#define GENERICTREENODE_H
 
#include "pup.h"
 
#include <map>
#include <vector>
#include <algorithm>
#include <set>
#include <list>
 
#include "OrientedBox.h"
#include "MultipoleMoments.h"
#include "keytype.h"
#include "GravityParticle.h"
 
namespace Tree {

  // C++11 syntax
  // using NodeKey = KeyType;
  typedef KeyType NodeKey;
  static const int NodeKeyBits = 8*sizeof(NodeKey);

  enum NodeType {
    Invalid = 1,
    Bucket,
    Internal,
    Boundary,
    NonLocal,
    Empty,
    Top,
    NonLocalBucket,
    Cached,
    CachedBucket,
    CachedEmpty
  };

  extern int numEmptyNodes;
 
class NodePool;

  class GenericTreeNode {
#ifdef CHANGA_REFACTOR_WALKCHECK
    public:
    bool touched;
    char by;
#endif
  protected:
    NodeType myType;
    NodeKey key;
 
    GenericTreeNode() : myType(Invalid), key(0), parent(0), firstParticle(0),
    lastParticle(0), remoteIndex(0), iParticleTypes(0), nSPH(0) {
#if COSMO_STATS > 0
      used = false;
#endif
#if INTERLIST_VER > 0
      numBucketsBeneath=0;
      startBucket=-1;
#ifdef CUDA
      nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
#endif
#ifdef CHANGA_REFACTOR_WALKCHECK
      touched = false;
      by = 'I';
#endif
    }
 
  public:
#if COSMO_STATS > 0
    bool used;
#endif

    MultipoleMoments moments;
    GenericTreeNode* parent;
    OrientedBox<cosmoType> boundingBox;
    OrientedBox<double> bndBoxBall;
    unsigned int iParticleTypes;
    int64_t nSPH;
    int firstParticle;
    int lastParticle;
    int remoteIndex;
    unsigned int particleCount;
    GravityParticle *particlePointer;

    int rungs;
 
#if INTERLIST_VER > 0
    int numBucketsBeneath;
    int startBucket;
#ifdef CUDA
    int nodeArrayIndex;
    int bucketArrayIndex;
#endif
#endif
    Vector3D<double> centerSm;
    double sizeSm;
    double fKeyMax;
    int iRank;

    GenericTreeNode(NodeKey k, NodeType type, int first, int last, GenericTreeNode *p) : myType(type), key(k), parent(p), firstParticle(first), lastParticle(last), remoteIndex(0) {
#if INTERLIST_VER > 0
      numBucketsBeneath=0;
      startBucket=-1;
#ifdef CUDA
      nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
#endif
    }
 
    virtual ~GenericTreeNode() { }
    virtual void fullyDelete() = 0;

    inline NodeType getType() const { return myType; }
    inline void setType(NodeType t) { myType = t; }

    inline NodeKey getKey() const { return key; }

    virtual unsigned int numChildren() const = 0;
    virtual GenericTreeNode* getChildren(int) = 0;
    virtual void setChildren(int, GenericTreeNode*) = 0;
    virtual NodeKey getChildKey(int) = 0;
    virtual NodeKey getParentKey() = 0;
    virtual int whichChild(NodeKey childkey) = 0;
#if INTERLIST_VER > 0
    virtual bool contains(NodeKey nodekey) = 0;
#endif

    bool isValid(){
      return (myType != Invalid);
    }

    bool isCached(){
      return (myType == Cached ||
              myType == CachedBucket ||
              myType == CachedEmpty);
    }

    bool isBucket(){
      return (myType == Bucket ||
              myType == CachedBucket ||
              myType == NonLocalBucket);
    }

    virtual void makeOctChildren(GravityParticle *part, int totalPart, int level, NodePool *pool = NULL) = 0;
    virtual void makeOrbChildren(GravityParticle *part, int totalPart, int level, int rootsLevel, bool (*compFnPtr[])(GravityParticle, GravityParticle), bool spatial, NodePool *pool = NULL) = 0;

    virtual void getChunks(int num, NodeKey *&ret) = 0;

    inline void makeBucket(GravityParticle *part) {
      myType = Bucket;
      iRank = CkMyRank();
#if INTERLIST_VER > 0
      numBucketsBeneath = 1;
#endif
      calculateRadiusBox(moments, boundingBox); /* set initial size */
      if(moments.getRadius() <= 0.0) {
          if(particleCount > 1)
              calculateRadiusFirstParticle(moments, &part[firstParticle],
                                           &part[lastParticle+1]);
          else
              moments.setRadius(1.0); // single particle boxes don't
                                      // need scaling.
      }
      boundingBox.reset();
      bndBoxBall.reset();
      iParticleTypes = 0;
      nSPH = 0;
      rungs = 0;
      for (int i = firstParticle; i <= lastParticle; ++i) {
        moments += part[i];
        boundingBox.grow(part[i].position);
    if(TYPETest(&part[i], TYPE_GAS)) {
        double fBallMax = part[i].fBallMax();
        bndBoxBall.grow(part[i].position
                + Vector3D<double>(fBallMax, fBallMax, fBallMax));
        bndBoxBall.grow(part[i].position
                - Vector3D<double>(fBallMax, fBallMax, fBallMax));
            nSPH++;
        }
    iParticleTypes |= part[i].iType;
        if (part[i].rung > rungs) rungs = part[i].rung;
      }
      if(particleCount > 1)
      calculateRadiusFarthestParticle(moments, &part[firstParticle],
                      &part[lastParticle+1]);
    }

    inline void makeEmpty() {
      myType = Empty;
      particleCount = 0;
#if INTERLIST_VER > 0
      numBucketsBeneath = 0;
#endif
      moments.clear();
      boundingBox.reset();
      bndBoxBall.reset();
      iParticleTypes = 0;
      nSPH = 0;
    }

    void getGraphViz(std::ostream &out);

    virtual NodeKey getLongestCommonPrefix(NodeKey k1, NodeKey k2)
    {
      CkAbort("getLongestCommonPrefix not implemented\n");
      return 0;
    }

    virtual int getLevel(NodeKey k) = 0;
    virtual GenericTreeNode *clone() const = 0;

    virtual void pup(PUP::er &p, int depth) = 0;
    virtual void pup(PUP::er &p) {
      int iType;
      if(p.isUnpacking()) {
        p | iType;
        myType = (NodeType) iType;
#ifdef CUDA
        // so that newly shipped nodes are not mistakenly assumed
        // to be present on the GPU
        nodeArrayIndex = -1;
      bucketArrayIndex = -1;
#endif
      } else {
        iType = (int) myType;
        p | iType;
      }
      p | key;
      p | moments;
      p | boundingBox;
      p | bndBoxBall;
      p | iParticleTypes;
      p | nSPH;
      p | firstParticle;
      p | lastParticle;
      p | remoteIndex;
      p | particleCount;
#if INTERLIST_VER > 0
      p | numBucketsBeneath;
      p | startBucket;
#endif
      p | centerSm;
      p | sizeSm;
      p | fKeyMax;
#ifdef CHANGA_REFACTOR_WALKCHECK
      p | touched;
      p | by;
#endif
    }
  };
 
class BinaryTreeNode;

class NodePool {
    int next;
    int szPool;
    std::list<BinaryTreeNode *> pools;
public:
    NodePool() {
        next = szPool = 128;  // Determines the size of the pools
    }
    ~NodePool();
    BinaryTreeNode *alloc_one();
    BinaryTreeNode *alloc_one(NodeKey k, NodeType type, int first,
                  int nextlast, BinaryTreeNode *p);
 
    };

  typedef std::map<NodeKey, GenericTreeNode *> NodeLookupType;

  class BinaryTreeNode : public GenericTreeNode {
  protected:
  public:
    BinaryTreeNode* children[2];
 
    BinaryTreeNode() : GenericTreeNode() {
      children[0] = 0;
      children[1] = 0;
    }
 
  public:
    BinaryTreeNode(NodeKey k, NodeType type, int first, int nextlast, BinaryTreeNode *p) : GenericTreeNode(k, type, first, nextlast, p) {
      children[0] = 0;
      children[1] = 0;
    }
 
    void fullyDelete() {
      if (children[0] != NULL) children[0]->fullyDelete();
      delete children[0];
      if (children[1] != NULL) children[1]->fullyDelete();
      delete children[1];
    }
 
    unsigned int numChildren() const {
      return 2;
    }
 
    GenericTreeNode* getChildren(int i) {
      CkAssert(i>=0 && i<2);
      return children[i];
    }
 
    void setChildren(int i, GenericTreeNode* node) {
      CkAssert(i>=0 && i<2);
      children[i] = (BinaryTreeNode*)node;
    }
 
    NodeKey getChildKey(int i) {
      CkAssert(i>=0 && i<2);
      return (key<<1) + i;
    }
 
    NodeKey getParentKey() {
      return (key>>1);
    }
 
    int getLevel(NodeKey k) {
      int i = 0;
      k >>= 1;
      while (k!=0) {
    k >>= 1;
    ++i;
      }
      return i;
    }
 
    NodeKey getLongestCommonPrefix(NodeKey k1, NodeKey k2){
 
      int l1 = getLevel(k1)+1;
      int l2 = getLevel(k2)+1;
 
      NodeKey a1 = k1 << (NodeKeyBits-l1);
      NodeKey a2 = k2 << (NodeKeyBits-l2);
 
      NodeKey a = a1 ^ a2;
      int i = 0;
      while( a < (NodeKey(1)<<(NodeKeyBits-1)) )
      {
        a = a << 1;
        i++;
      }
      return (k1 >> (l1-i)); // or k2 >> (l2-i)
    };
 
 
    int whichChild(NodeKey child) {
      int thisLevel = getLevel(key);
      int childLevel = getLevel(child);
      return ((child >> (childLevel-thisLevel-1)) ^ (key << 1));
    }
 
#if INTERLIST_VER > 0
    bool contains(NodeKey node) {
      int thisLevel = getLevel(key);
      int nodeLevel = getLevel(node);
 
      if(nodeLevel<thisLevel) return false;
 
      return ((node>>(nodeLevel-thisLevel))==key);
    }
#endif
 
    bool isLeftChild() const {
      return (dynamic_cast<BinaryTreeNode *>(parent) && dynamic_cast<BinaryTreeNode *>(parent)->children[0] == this);
    }
 
    bool isRightChild() const {
      return (dynamic_cast<BinaryTreeNode *>(parent) && dynamic_cast<BinaryTreeNode *>(parent)->children[1] == this);
    }
 
    BinaryTreeNode* getSibling() const {
      BinaryTreeNode* p = dynamic_cast<BinaryTreeNode *>(parent);
      if(p)
    return (p->children[0] == this ? p->children[1] : p->children[0]);
      else
    return 0;
    }

    void makeOctChildren(GravityParticle *part, int totalPart, int level,
             NodePool *pool=NULL) {
      if(pool) {
      children[0] = pool->alloc_one();
      children[1] = pool->alloc_one();
      }
      else {
      children[0] = new BinaryTreeNode();
      children[1] = new BinaryTreeNode();
      }
      children[0]->parent = this;
      children[1]->parent = this;
      children[0]->boundingBox = boundingBox;
      children[1]->boundingBox = boundingBox;
      double split;
      switch (level%3) {
      case 0:
    split = 0.5 * (boundingBox.greater_corner.x + boundingBox.lesser_corner.x);
    children[0]->boundingBox.greater_corner.x = split;
    children[1]->boundingBox.lesser_corner.x = split;
    break;
      case 1:
    split = 0.5 * (boundingBox.greater_corner.y + boundingBox.lesser_corner.y);
    children[0]->boundingBox.greater_corner.y = split;
    children[1]->boundingBox.lesser_corner.y = split;
    break;
      case 2:
    split = 0.5 * (boundingBox.greater_corner.z + boundingBox.lesser_corner.z);
    children[0]->boundingBox.greater_corner.z = split;
    children[1]->boundingBox.lesser_corner.z = split;
    break;
      }
      children[0]->key = key << 1;
      children[1]->key = (key << 1) + 1;
      children[0]->firstParticle = firstParticle;
      children[1]->lastParticle = lastParticle;
 
      // mask holds the bit that determines a particle's membership in
      // either the left (0) or right (1) child.
      SFC::Key mask = SFC::Key(1) << ((SFC::KeyBits-1) - level);
      SFC::Key leftBit = part[firstParticle].key & mask;
      SFC::Key rightBit = part[lastParticle].key & mask;
 
      if (leftBit == rightBit) {
    // we have only one child, the other is either NonLocal or Empty
    if (leftBit) {
      // left child missing
      if (firstParticle == 0) {
              // We are at the left domain boundary: the left child is
              // completely on another TreePiece
              children[0]->myType = NonLocal;
              children[1]->myType = Boundary;
      } else {
        children[0]->makeEmpty();
        children[1]->myType = lastParticle==totalPart+1 ? Boundary : Internal;
      }
      children[0]->lastParticle = firstParticle-1;
          children[0]->particleCount = 0;
      children[1]->firstParticle = firstParticle;
      children[1]->particleCount = particleCount;
    } else {
      // right child missing
      if (lastParticle == totalPart+1) {
              // We are at the right domain boundary: the right child is
              // completely on another TreePiece
              children[1]->myType = NonLocal;
              children[0]->myType = Boundary;
      } else {
        children[1]->makeEmpty();
        children[0]->myType = firstParticle==0 ? Boundary : Internal;
      }
      children[1]->firstParticle = lastParticle+1;
          children[1]->particleCount = 0;
      children[0]->lastParticle = lastParticle;
      children[0]->particleCount = particleCount;
    }
      } else if (leftBit < rightBit) {
    // both children are present
    if (firstParticle == 0 && leftBit != (part[1].key & mask)) {
          // the left child is NonLocal: the boundary particle (from
          // another TP) is in the left child, but the rest of the
          // particles are in the right child.
      children[0]->myType = NonLocal;
      children[0]->lastParticle = firstParticle;
      children[1]->myType = lastParticle==totalPart+1 ? Boundary : Internal;
      children[1]->firstParticle = firstParticle+1;
      children[1]->particleCount = particleCount;
    } else if (lastParticle == totalPart+1 && rightBit != (part[totalPart].key & mask)) {
          // the right child is NonLocal: the boundary particle (from
          // another TP) is in the right child, but the rest of the
          // particles are in the left child.
      children[1]->myType = NonLocal;
      children[1]->firstParticle = lastParticle;
      children[0]->myType = firstParticle==0 ? Boundary : Internal;
      children[0]->lastParticle = lastParticle-1;
      children[0]->particleCount = particleCount;
    } else {
          // the splitting point is somewhere in the middle
          // The following statement finds the first particle with
          // splitting bit (see the variable 'mask' above) set.
      GravityParticle *splitParticle = std::lower_bound(&part[firstParticle],
            &part[lastParticle+1],
            (GravityParticle)(part[lastParticle].key
                      & ((~ (SFC::Key)0) << ((SFC::KeyBits-1)-level))));
      children[0]->lastParticle = splitParticle - part - 1;
      children[1]->firstParticle = splitParticle - part;
      children[0]->particleCount = children[0]->lastParticle - firstParticle + 1;
      children[1]->particleCount = lastParticle - children[1]->firstParticle + 1;
      children[0]->myType = children[0]->particleCount==0 ? Empty : (firstParticle==0 ? Boundary : Internal);
      children[1]->myType = children[1]->particleCount==0 ? Empty : (lastParticle==totalPart+1 ? Boundary : Internal);
          if (children[0]->myType==Empty) children[0]->makeEmpty();
          if (children[1]->myType==Empty) children[1]->makeEmpty();
      if (firstParticle == 0) children[0]->particleCount --;
      if (lastParticle == totalPart+1) children[1]->particleCount --;
    }
      } else {
    CkAbort("Ok, the particles should be ordered so, how the hell did we get here?!?");
      }
      children[0]->particlePointer = &part[children[0]->firstParticle];
      children[1]->particlePointer = &part[children[1]->firstParticle];
#if INTERLIST_VER > 0
      if(children[0]->myType == NonLocal) children[0]->numBucketsBeneath = 0;
      if(children[1]->myType == NonLocal) children[1]->numBucketsBeneath = 0;
      // empty nodes are makeEmpty()'ed, so that the numbucketsbeneath them are 0
#endif
    }
 
    // Constructs 2 children of this node based on ORB decomposition
    void makeOrbChildren(GravityParticle *part, int totalPart, int level,
                 int rootsLevel,
             bool (*compFnPtr[])(GravityParticle, GravityParticle),
             bool spatial, NodePool *pool=NULL) {
 
      if(pool) {
      children[0] = pool->alloc_one();
      children[1] = pool->alloc_one();
      }
      else {
          children[0] = new BinaryTreeNode();
          children[1] = new BinaryTreeNode();
      }
      children[0]->parent = this;
      children[1]->parent = this;
      children[0]->boundingBox = boundingBox;
      children[1]->boundingBox = boundingBox;
 
      children[0]->key = key << 1;
      children[1]->key = (key << 1) + 1;
      children[0]->firstParticle = firstParticle;
      children[1]->lastParticle = lastParticle;
 
      //CkPrintf("children keys:%lld,%lld\n",children[0]->key,children[1]->key);
      if(level<rootsLevel){
        //This branch is taken for levels above the TreePiece root level
        //TreePiece root level is the level from where onwards data is internal
        SFC::Key tmp = 1 << (level+1);
        tmp = children[0]->key - tmp;
        SFC::Key tmp2 = part[0].key >> (rootsLevel-level-1);
 
        if(level+1 != rootsLevel){
          if(tmp==tmp2){
            children[0]->myType = Boundary;
            children[1]->myType = NonLocal;
              children[1]->firstParticle = lastParticle+1;
              children[0]->lastParticle = lastParticle;
              children[0]->particleCount = particleCount;
          }
          else{
            children[0]->myType = NonLocal;
            children[1]->myType = Boundary;
              children[0]->lastParticle = firstParticle-1;
              children[1]->firstParticle = firstParticle;
              children[1]->particleCount = particleCount;
          }
        }
        else{ //Special case for TreePiece root level
          if(tmp==tmp2){
            children[0]->myType = Internal;
            children[1]->myType = NonLocal;
              children[1]->firstParticle = lastParticle+1;
              children[0]->lastParticle = lastParticle;
              children[0]->particleCount = particleCount;
              if(firstParticle==0)
              children[0]->firstParticle = firstParticle+1;
              if(lastParticle==totalPart+1)
              children[0]->lastParticle = lastParticle-1;
          }
          else{
            children[0]->myType = NonLocal;
            children[1]->myType = Internal;
              children[0]->lastParticle = firstParticle-1;
              children[1]->firstParticle = firstParticle;
              children[1]->particleCount = particleCount;
              if(firstParticle==0)
              children[1]->firstParticle = firstParticle+1;
              if(lastParticle==totalPart+1)
              children[1]->lastParticle = lastParticle-1;
          }
        }
      }
      else{ //Below the TreePiece root level
        float len=0.0,len2=0.0;
        int dim;
    
    if(spatial) { // "squeeze" box before division
        boundingBox.reset();
        for (int i = firstParticle; i <= lastParticle; ++i)
        boundingBox.grow(part[i].position);
        }
 
        len=boundingBox.greater_corner.x-boundingBox.lesser_corner.x;
        dim=0;
        CkAssert(len>=0.0);
        len2=boundingBox.greater_corner.y-boundingBox.lesser_corner.y;
        CkAssert(len2>=0.0);
        if(len2>len) { len = len2; dim=1; }
        len2=boundingBox.greater_corner.z-boundingBox.lesser_corner.z;
        CkAssert(len2>=0.0);
        if(len2>len) { len = len2; dim=2; }
 
        //Sort the particles in longest dimension
        std::sort(&part[firstParticle],&part[lastParticle+1],compFnPtr[dim]);
 
        //Middle particle is the splitting point to get 2 children
 
        //Compress the bounding box of this node
        boundingBox.greater_corner[dim] = part[lastParticle].position[dim];
        boundingBox.lesser_corner[dim] = part[firstParticle].position[dim];
 
    if(!spatial) {
        if((lastParticle-firstParticle+1)%2==0){
          children[0]->lastParticle = firstParticle + (lastParticle-firstParticle-1)/2;
          children[1]->firstParticle = firstParticle + (lastParticle-firstParticle+1)/2;
          children[0]->particleCount = (lastParticle-firstParticle+1)/2;
          children[1]->particleCount = (lastParticle-firstParticle+1)/2;
        }
        else{
          children[0]->lastParticle = firstParticle + (lastParticle-firstParticle)/2;
          children[1]->firstParticle = firstParticle + (lastParticle-firstParticle+2)/2;
          children[0]->particleCount = (lastParticle-firstParticle+2)/2;
          children[1]->particleCount = (lastParticle-firstParticle)/2;
        }
    }
    else {
        GravityParticle dummy;  // dummy for splitting
        GravityParticle* divStart = &part[firstParticle];
        Vector3D<double> divide(0.0,0.0,0.0);
        divide[dim] = part[firstParticle].position[dim] + 0.5*len;
        dummy.position = divide;
        GravityParticle* divEnd
        = std::upper_bound(&part[firstParticle],&part[lastParticle+1],
                   dummy,compFnPtr[dim]);
        children[0]->lastParticle = firstParticle + (divEnd - divStart) - 1;
        children[1]->firstParticle = firstParticle + (divEnd - divStart);
        children[0]->particleCount = divEnd - divStart;
        children[1]->particleCount = 1 + (lastParticle-firstParticle)
        - (divEnd - divStart);
        CkAssert(children[0]->particleCount > 0);
        CkAssert(children[1]->particleCount > 0);
    }
 
        children[0]->myType = Internal;
        children[1]->myType = Internal;
 
        //Compress the bounding boxes of both children too
        children[0]->boundingBox.lesser_corner = boundingBox.lesser_corner;
        children[0]->boundingBox.greater_corner = boundingBox.greater_corner;
        children[0]->boundingBox.greater_corner[dim] = part[children[0]->lastParticle].position[dim];
 
        children[1]->boundingBox.lesser_corner = boundingBox.lesser_corner;
        children[1]->boundingBox.lesser_corner[dim] = part[children[1]->firstParticle].position[dim];
        children[1]->boundingBox.greater_corner = boundingBox.greater_corner;
 
      }
      children[0]->particlePointer = &part[children[0]->firstParticle];
      children[1]->particlePointer = &part[children[1]->firstParticle];
 
    }
 
    // implemented in the .cpp
    GenericTreeNode *clone() const;

    void getChunks(int num, NodeKey *&ret) {
      int i = 0;
      int base = num;
      while (base > 1) {
    base >>= 1;
    i++;
      }
      base = 1 << i;
      // base now contains the greatest power of two less or equal to num
      int additional = num - base;
      if (ret==NULL) ret = new NodeKey[num];
      for (int j=additional, k=additional*2; j<base; ++j, ++k) ret[k] = j + base;
      base <<= 1;
      for (int j=0; j<additional*2; ++j) ret[j] = j + base;
 
    }
 
    int countDepth(int depth) {
      int count = 1;
      if (depth != 0) {
        if (children[0] != NULL) count += children[0]->countDepth(depth-1);
        if (children[1] != NULL) count += children[1]->countDepth(depth-1);
      }
      return count;
    }
 
    // extraSpace is used by the CacheInterface to store a pointer to
    // the message so it can be freed when we are finished.
    int packNodes(BinaryTreeNode *buffer, int depth, int extraSpace=0) {
      //CkPrintf("Entering packNodes: this=%p, buffer=%p, depth=%d\n",this,buffer,depth);
      //*buffer = *this;
      memcpy((void *)buffer, (void *)this, sizeof(*this));
      buffer->parent = NULL;
      buffer->particlePointer = NULL;
#if INTERLIST_VER > 0 && defined CUDA
      buffer->nodeArrayIndex = -1;
      buffer->bucketArrayIndex = -1;
#endif
      int used = 1;
      if (depth != 0) {
        if (children[0] != NULL) {
          BinaryTreeNode *nextBuf = (BinaryTreeNode *) (((char*)buffer) + used * ALIGN_DEFAULT(sizeof(BinaryTreeNode)+extraSpace));
          buffer->children[0] = (BinaryTreeNode*)(((char*)nextBuf) - ((char*)buffer));
          //CkPrintf("Entering child 0: offset %ld\n",buffer->children[0]);
          used += children[0]->packNodes(nextBuf, depth-1, extraSpace);
        } else {
          //CkPrintf("Excluding child 0\n");
          buffer->children[0] = NULL;
        }
        if (children[1] != NULL) {
          BinaryTreeNode *nextBuf = (BinaryTreeNode *) (((char*)buffer) + used * ALIGN_DEFAULT(sizeof(BinaryTreeNode)+extraSpace));
          buffer->children[1] = (BinaryTreeNode*)(((char*)nextBuf) - ((char*)buffer));
          //CkPrintf("Entering child 1: offset %ld\n",buffer->children[1]);
          used += children[1]->packNodes(nextBuf, depth-1, extraSpace);
        } else {
          //CkPrintf("Excluding child 1\n");
          buffer->children[1] = NULL;
        }
      } else {
        //CkPrintf("Depth reached\n");
        buffer->children[0] = buffer->children[1] = NULL;
      }
      //CkAssert((long int)buffer->children[0] < 10000);
      //CkAssert((long int)buffer->children[1] < 10000);
      //CkPrintf("Returning used = %d\n",used);
      return used;
    }
 
    void unpackNodes() {
      if (children[0] != NULL) {
        //CkAssert((long int)children[0] < 10000);
        children[0] = (BinaryTreeNode*)(((long int)children[0]) + ((char*)this));
        children[0]->parent = this;
        children[0]->unpackNodes();
      }
      if (children[1] != NULL) {
        //CkAssert((long int)children[1] < 10000);
        children[1] = (BinaryTreeNode*)(((long int)children[1]) + ((char*)this));
        children[1]->parent = this;
        children[1]->unpackNodes();
      }
    }
 
    void pup(PUP::er &p) { pup(p, -1); }
    void pup(PUP::er &p, int depth);/* {
      //CkPrintf("Pupper of BinaryTreeNode(%d) called for %s (%d)\n",depth,p.isPacking()?"Packing":p.isUnpacking()?"Unpacking":"Sizing",p.isSizing()?((PUP::sizer*)&p)->size():((PUP::mem*)&p)->size());
      GenericTreeNode::pup(p);
      int isNull;
      for (int i=0; i<2; ++i) {
    isNull = (children[i]==NULL || depth==0) ? 0 : 1;
    p | isNull;
    CkAssert(isNull==0 || isNull==1);
    if (isNull != 0 && depth != 0) {
      if (p.isUnpacking()) children[i] = new BinaryTreeNode();
      children[i]->pup(p, depth-1);
      if (p.isUnpacking()) children[i]->parent = this;
    }
      }
      }*/
  };
 
inline
NodePool::~NodePool() {
        while(!pools.empty()) {
            delete [] pools.front();
            pools.pop_front();
            }
        }
 
inline BinaryTreeNode *
NodePool::alloc_one() {
        if(next >= szPool) { // need new pool block
            BinaryTreeNode *nodes = new BinaryTreeNode[szPool];
            next = 0;
            pools.push_front(nodes);
            }
        return &pools.front()[next++];
        }
 
inline BinaryTreeNode *
NodePool::alloc_one(NodeKey k, NodeType type, int first, int nextlast,
            BinaryTreeNode *p) {
    BinaryTreeNode *one = alloc_one();
    return new (one) BinaryTreeNode(k, type, first, nextlast, p);
    }

  class OctTreeNode : public GenericTreeNode {
  protected:
  public:
    OctTreeNode* children[8];
 
    OctTreeNode() : GenericTreeNode() {
      for (int i = 0; i < 8; i++) {
        children[i] = 0;
      }
    }
 
  public:
    OctTreeNode(NodeKey k, NodeType type, int first, int last, BinaryTreeNode *p) : GenericTreeNode(k, type, first, last, p) {
      for (int i = 0; i < 8; i++) {
        children[i] = 0;
      }
    }
 
    void fullyDelete() {
      for (int i = 0; i < numChildren(); i++) {
        if (children[i] != NULL) {
          children[i]->fullyDelete();
          delete children[i];
        }
      }
    }
 
    virtual unsigned int numChildren() const {
      return 8;
    }
 
    virtual GenericTreeNode* getChildren(int i) {
      CkAssert(i>=0 && i<8);
      return children[i];
    }
 
    virtual void setChildren(int i, GenericTreeNode* node) {
      CkAssert(i>=0 && i<8);
      children[i] = (OctTreeNode*)node;
    }
 
    NodeKey getChildKey(int i) {
      CkAssert(i>=0 && i<8);
      return (key<<3) + i;
    }
 
    NodeKey getParentKey() {
      return (key>>3);
    }
 
    int whichChild(NodeKey child) {
      return (child ^ (key<<3));
    }
 
    int getLevel(NodeKey k) {
      int i = 0;
      k >>= 3;
      while (k!=0) {
    k >>= 3;
    ++i;
      }
      return i;
    }
 
#if INTERLIST_VER > 0
    bool contains(NodeKey node) {
 
      return true;
    }
#endif
 
    void makeOctChildren(GravityParticle *part, int totalPart, int level,
             NodePool *pool = NULL) {}
 
    void makeOrbChildren(GravityParticle *part, int totalPart, int level, int rootsLevel, bool (*compFnPtr[])(GravityParticle, GravityParticle), bool spatial,
             NodePool *pool = NULL) {}
 
    GenericTreeNode *clone() const {
      OctTreeNode *tmp = new OctTreeNode();
      *tmp = *this;
      return tmp;
    }
 
    /*
    int getNumChunks(int num) {
      int i = 0;
      while (num > 1) {
    num >>= 3;
    i++;
      }
      return 1 << (3*i);
    }
    */
 
    void getChunks(int num, NodeKey *&ret) {
      return;
    }
 
    void pup(PUP::er &p);
    void pup(PUP::er &p, int depth);
  };

  enum GenericTrees {
    Binary_Oct,
    Oct_Oct,
    Binary_ORB
  };

 
  class compare{ //Defines the comparison operator on the map used in balancer
  public:
    compare(){}
 
    bool operator()(NodeKey key1, NodeKey key2) const {
 
      NodeKey tmp = NodeKey(1);
      int len1=0, len2=0;
      int cnt=1;
 
      while(tmp<=key1 || tmp<=key2){
        tmp<<=1;
        cnt++;
        if(len1==0 && tmp>key1){
          len1=cnt-1;
        }
        if(len2==0 && tmp>key2){
          len2=cnt-1;
        }
      }
 
      if(len1==len2){
        return key1<key2;
      }
      else if(len1>len2){
        key1>>=(len1-len2);
        return key1<key2;
      }
      else{
        key2>>=(len2-len1);
        return key1<key2;
      }
    }
  };

  inline void operator|(PUP::er &p, NodeType &nt) {
    int nti;
    if (p.isUnpacking()) {
      p | nti;
      nt = (NodeType)nti;
    } else {
      nti = (int)nt;
      p | nti;
    }
  }

  inline void operator|(PUP::er &p, GenericTrees &gt) {
    int gti;
    if (p.isUnpacking()) {
      p | gti;
      gt = (GenericTrees)gti;
    } else {
      gti = (int)gt;
      p | gti;
    }
  }
 
} //close namespace Tree
 
#endif //GENERICTREENODE_H