CarpetLib: New class bboxset2

Rename bboxset to bboxset1.

Implement new class bboxset2, which uses a different internal datastructure than bboxset1.

Make bboxset a typedef for either bboxset1 (default) or bboxset2, as selected by the compile-time macro CARPET_BBOXSET2.

Disable all bboxset2 code if CARPET_NO_BBOXSET2 is given, since bboxset2 uses newer C++ constructs not available on some older compilers.

10 files changed, 2786 insertions, 1226 deletions

diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index e94f17f43..e4970669a 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -1,794 +1 @@
-#include <algorithm>
-#include <cassert>
-#include <iostream>
-#include <limits>
-#include <set>
-#include <stack>
-
-#include "defs.hh"
-
 #include "bboxset.hh"
-
-using namespace std;
-
-
-
-#define SKIP_NORMALIZE(s) skip_normalize_t skip_normalize_v(s)
-// #define SKIP_NORMALIZE(s) do { } while(0)
-
-
-
-// Constructors
-template<typename T, int D>
-bboxset<T,D>::bboxset ()
-  : skip_normalize(false)
-{
-  assert (invariant());
-}
-
-template<typename T, int D>
-bboxset<T,D>::bboxset (const box& b)
-  : skip_normalize(false)
-{
-  //S if (not b.empty()) bs.insert(b);
-  if (not b.empty()) bs.push_back(b);
-  assert (invariant());
-}
-
-template<typename T, int D>
-bboxset<T,D>::bboxset (const bboxset& s)
-  : bs(s.bs), skip_normalize(false)
-{
-  assert (invariant());
-}
-
-template<typename T, int D>
-bboxset<T,D>::bboxset (const list<box>& lb)
-  : skip_normalize(false)
-{
-  SKIP_NORMALIZE(*this);
-  for (typename list<box>::const_iterator
-         li = lb.begin(), le = lb.end(); li != le; ++ li)
-  {
-    *this |= *li;
-  }
-}
-
-template<typename T, int D>
-bboxset<T,D>::bboxset (const vector<box>& lb)
-  : skip_normalize(false)
-{
-  SKIP_NORMALIZE(*this);
-  for (typename vector<box>::const_iterator
-         vi = lb.begin(), ve = lb.end(); vi != ve; ++ vi)
-  {
-    *this |= *vi;
-  }
-}
-
-template<typename T, int D>
-bboxset<T,D>::bboxset (const vector<list<box> >& vlb)
-  : skip_normalize(false)
-{
-  SKIP_NORMALIZE(*this);
-  for (typename vector<list<box> >::const_iterator
-         vli = vlb.begin(), vle = vlb.end(); vli != vle; ++ vli)
-  {
-    *this |= bboxset (*vli);
-  }
-}
-
-template<typename T, int D>
-template<typename U>
-bboxset<T,D>::bboxset (const vector<U>& vb, const bbox<T,D> U::* const v)
-  : skip_normalize(false)
-{
-  SKIP_NORMALIZE(*this);
-  for (typename vector<U>::const_iterator
-         vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi)
-  {
-    *this |= (*vi).*v;
-  }
-}
-
-template<typename T, int D>
-template<typename U>
-bboxset<T,D>::bboxset (const vector<U>& vb, const bboxset U::* const v)
-  : skip_normalize(false)
-{
-  SKIP_NORMALIZE(*this);
-  for (typename vector<U>::const_iterator
-         vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi)
-  {
-    *this |= (*vi).*v;
-  }
-}
-
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::poison () {
-  return bboxset (bbox<T,D>::poison());
-}
-
-
-
-// Invariant
-template<typename T, int D>
-bool bboxset<T,D>::invariant () const {
-  // This is very slow when there are many bboxes
-#if 0 && defined(CARPET_DEBUG)
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    if ((*bi).empty()) return false;
-    if (not (*bi).is_aligned_with(*bs.begin())) return false;
-    // check for overlap (quadratic -- expensive)
-    for (const_iterator bi2=begin(); bi2!=bi; ++bi2) {
-      if ((*bi2).intersects(*bi)) return false;
-    }
-  }
-#endif
-  return true;
-}
-
-
-
-// Normalisation
-template<typename T, int D>
-void bboxset<T,D>::normalize ()
-{
-  if (skip_normalize) return;
-  assert (invariant());
-  
-  bboxset const oldbs = * this;
-  size_type const oldsize = this->size();
-  
-  // Split all bboxes into small pieces which have all their
-  // boundaries aligned.
-  for (int d=0; d<D; ++d) {
-    // Find all boundaries
-    //S typedef set<T> buf;
-    typedef vector<T> buf;
-    buf sbnds;
-    sbnds.reserve (2 * bs.size());
-    for (typename bset::const_iterator
-           si = bs.begin(), se = bs.end(); si != se; ++ si)
-    {
-      box const & b = * si;
-      int const bstr = b.stride()[d];
-      int const blo = b.lower()[d];
-      int const bhi = b.upper()[d] + bstr;
-      //S sbnds.insert (blo);
-      //S sbnds.insert (bhi);
-      sbnds.push_back (blo);
-      sbnds.push_back (bhi);
-    }
-    sort (sbnds.begin(), sbnds.end());
-    typename buf::iterator const last = unique (sbnds.begin(), sbnds.end());
-    sbnds.resize (last - sbnds.begin());
-    // Split bboxes
-    bset nbs;
-    for (typename bset::const_iterator
-           si = bs.begin(), se = bs.end(); si != se; ++ si)
-    {
-      box const & b = * si;
-      int const bstr = b.stride()[d];
-      int const blo = b.lower()[d];
-      int const bhi = b.upper()[d] + bstr;
-      typename buf::const_iterator const ilo
-        = lower_bound (sbnds.begin(), sbnds.end(), blo);
-      typename buf::const_iterator const ihi
-        = lower_bound (sbnds.begin(), sbnds.end(), bhi);
-      assert (ilo != sbnds.end());
-      assert (ihi != sbnds.end());
-      assert (* ilo == blo);
-      assert (* ihi == bhi);
-      // Split one bbox
-      for (typename buf::const_iterator curr = ilo; curr != ihi; ++ curr) {
-        typename buf::const_iterator next = curr;
-        advance (next, 1);
-        int const nblo = * curr;
-        int const nbhi = * next;
-        assert (nbhi > nblo);   // ensure that the set remains sorted
-        box const nb (b.lower().replace(d, nblo),
-                      b.upper().replace(d, nbhi - bstr),
-                      b.stride());
-        //S nbs.insert (nb);
-        nbs.push_back (nb);
-      }
-    }
-    // Replace old set
-    bs.clear ();
-    bs.swap (nbs);
-    assert (invariant());
-  }
-  
-  // Combine bboxes if possible
-  for (int d=0; d<D; ++d) {
-    bset nbs;
-    while (not bs.empty()) {
-      typename bset::iterator si = bs.begin();
-      assert (si != bs.end());
-      
-      box const b = * si;
-      int const bstr = b.stride()[d];
-      int const blo = b.lower()[d];
-      int const bhi = b.upper()[d] + bstr;
-      
-      for (typename bset::iterator
-             nsi = nbs.begin(), nse = nbs.end(); nsi != nse; ++ nsi)
-      {
-        box const nb = * nsi;
-        int const nblo = nb.lower()[d];
-        int const nbhi = nb.upper()[d] + bstr;
-        
-        box const mb (nb.lower().replace(d, blo),
-                      nb.upper().replace(d, bhi - bstr),
-                      nb.stride());
-        
-        // Check whether the other dimensions match
-        if (b == mb) {
-          // Check whether the bboxes are adjacent in this dimension
-          if (nbhi == blo) {
-            // Combine boxes, nb < b
-            box const cb (b.lower().replace(d, nblo),
-                          b.upper(),
-                          b.stride());
-            bs.erase (si);
-            nbs.erase (nsi);
-            //S bs.insert (cb);
-            bs.push_back (cb);
-            goto done;
-          } else if (bhi == nblo) {
-            // Combine boxes, b < nb
-            box const cb (b.lower(),
-                          b.upper().replace(d, nbhi - bstr),
-                          b.stride());
-            bs.erase (si);
-            nbs.erase (nsi);
-            //S bs.insert (cb);
-            bs.push_back (cb);
-            goto done;
-          }
-        }
-      }
-      bs.erase (si);
-      //S nbs.insert (b);
-      nbs.push_back (b);
-    done:;
-    }
-    bs.swap (nbs);
-    assert (invariant());
-  }
-  
-  size_type const newsize = this->size();
-  
-  // Can't use operators on *this since these would call normalize again
-  // assert (*this == oldbs);
-  assert (newsize == oldsize);
-}
-
-
-
-// Accessors
-// cost: O(n)
-template<typename T, int D>
-typename bboxset<T,D>::size_type bboxset<T,D>::size () const
-{
-  size_type s = 0;
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    const size_type bsz = (*bi).size();
-    assert (numeric_limits<size_type>::max() - bsz >= s);
-    s += bsz;
-  }
-  return s;
-}
-
-
-
-// Queries
-
-// Containment
-// cost: O(n)
-template<typename T, int D>
-bool bboxset<T,D>::contains (const vect<T,D>& x) const
-{
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    if ((*bi).contains(x)) return true;
-  }
-  return false;
-}
-
-// Intersection
-// cost: O(n)
-template<typename T, int D>
-bool bboxset<T,D>::intersects (const box& b) const
-{
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    if ((*bi).intersects(b)) return true;
-  }
-  return false;
-}
-
-
-
-// Add (bboxes that don't overlap)
-// cost: O(1)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator+= (const box& b)
-{
-  if (b.empty()) return *this;
-  // This is very slow when there are many bboxes
-#if 0 && defined(CARPET_DEBUG)
-  // check for overlap
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    assert (not (*bi).intersects(b));
-  }
-#endif
-  //S bs.insert(b);
-  bs.push_back(b);
-  assert (invariant());
-  normalize();
-  return *this;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s)
-{
-  SKIP_NORMALIZE(*this);
-  for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
-    *this += *bi;
-  }
-  assert (invariant());
-  return *this;
-}
-
-//L cost: O(1)
-// cost O(n)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::add_transfer (bboxset& s)
-{
-#ifdef BBOXSET_LIST
-  bs.splice (bs.end(), s.bs);
-#else
-  bs.insert (bs.end(), s.bs.begin(), s.bs.end());
-  s.bs.clear ();
-#endif
-  assert (invariant());
-  normalize();
-  return *this;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator+ (const box& b) const
-{
-  bboxset r(*this);
-  r += b;
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator+ (const bboxset& s) const
-{
-  bboxset r(*this);
-  r += s;
-  assert (r.invariant());
-  return r;
-}
-
-
-
-// Union
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator|= (const box& b)
-{
-  if (b.empty()) return *this;
-#if 0
-  // this has a cost of O(n^2)
-  bboxset tmp = b - *this;
-  add_transfer (tmp);
-#else
-  // this has a cost of O(n)
-  bset oldbs;
-  oldbs.swap(bs);
-  bs.push_back(b);
-  SKIP_NORMALIZE(*this);
-  for (const_iterator bi=oldbs.begin(), be=oldbs.end(); bi!=be; ++bi) {
-    bboxset tmp = *bi - b;
-    add_transfer(tmp);
-  }
-#endif
-  assert (invariant());
-  return *this;
-}
-
-// cost: O(n^2)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s)
-{
-  bboxset tmp = s - *this;
-  add_transfer (tmp);
-  assert (invariant());
-  return *this;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator| (const box& b) const
-{
-  bboxset r(*this);
-  r |= b;
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n^2)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator| (const bboxset& s) const
-{
-  bboxset r(*this);
-  r |= s;
-  assert (r.invariant());
-  return r;
-}
-
-
-
-// Intersection
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator& (const box& b) const
-{
-  // start with an empty set
-  bboxset r;
-  {
-    SKIP_NORMALIZE(r);
-    // walk all my elements
-    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-      // insert the intersection with the bbox
-      r += *bi & b;
-    }
-  }
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const
-{
-  // start with an empty set
-  bboxset r;
-  {
-    SKIP_NORMALIZE(r);
-    // walk all the bboxes
-    for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
-      // insert the intersection with this bbox
-      bboxset tmp = *this & *bi;
-      r.add_transfer (tmp);
-    }
-  }
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator&= (const box& b)
-{
-  *this = *this & b;
-  assert (invariant());
-  return *this;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s)
-{
-  *this = *this & s;
-  assert (invariant());
-  return *this;
-}
-
-
-
-// Difference
-// cost: O(1)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2)
-{
-  assert (b1.is_aligned_with(b2));
-  if (b1.empty()) return bboxset<T,D>();
-  if (not b1.intersects(b2)) return bboxset<T,D>(b1);
-  vect<T,D> const b2lo = max (b1.lower(), b2.lower());
-  vect<T,D> const b2up = min (b1.upper(), b2.upper());
-  vect<T,D> const & b1lo = b1.lower();
-  vect<T,D> const & b1up = b1.upper();
-  vect<T,D> const & str = b1.stride();
-  bboxset<T,D> r;
-  SKIP_NORMALIZE(r);
-  {
-    for (int d=0; d<D; ++d) {
-      // make resulting bboxes as large as possible in x-direction
-      // (for better consumption by Fortranly ordered arrays)
-      vect<T,D> lb, ub;
-      bbox<T,D> b;
-      for (int dd=0; dd<D; ++dd) {
-        if (dd<d) {
-          lb[dd] = b2lo[dd];
-          ub[dd] = b2up[dd];
-        } else if (dd>d) {
-          lb[dd] = b1lo[dd];
-          ub[dd] = b1up[dd];
-        }
-      }
-      lb[d] = b1lo[d];
-      ub[d] = b2lo[d] - str[d];
-      b = bbox<T,D>(lb,ub,str);
-      r += b;
-      lb[d] = b2up[d] + str[d];
-      ub[d] = b1up[d];
-      b = bbox<T,D>(lb,ub,str);
-      r += b;
-    }
-  }
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator- (const box& b) const
-{
-  // start with an empty set
-  bboxset r;
-  {
-    SKIP_NORMALIZE(r);
-    // walk all my elements
-    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-      // insert the difference with the bbox
-      bboxset tmp = *bi - b;
-      r.add_transfer (tmp);
-    }
-  }
-  assert (r.invariant());
-  return r;
-}
-  
-// cost: O(n^2)
-template<typename T, int D>
-bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s)
-{
-  SKIP_NORMALIZE(*this);
-  for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
-    *this -= *bi;
-  }
-  assert (invariant());
-  return *this;
-}
-
-// cost: O(n^2)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const
-{
-  bboxset r(*this);
-  r -= s;
-  assert (r.invariant());
-  return r;
-}
-
-// cost: O(n^2)
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  bboxset<T,D> r = bboxset<T,D>(b) - s;
-  assert (r.invariant());
-  return r;
-}
-
-
-
-// cost: O(n)
-template<typename T, int D>
-typename bboxset<T,D>::box bboxset<T,D>::container () const
-{
-  box b;
-  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-    b = b.expanded_containing(*bi);
-  }
-  return b;
-}
-
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::pseudo_inverse (const int n) const
-{
-  assert (not empty());
-  box const c = container().expand(n,n);
-  return c - *this;
-}
-
-// cost: O(n^2) in general, but only O(n) for shifting
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi)
-  const
-{
-  bboxset res;
-  {
-    SKIP_NORMALIZE(res);
-    if (all (lo == -hi)) {
-      // Special case for shifting, since this is faster
-      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-        res += (*bi).expand(lo,hi);
-      }
-    } else {
-      // We don't know (yet?) how to shrink a set
-      assert (all (lo>=0 and hi>=0));
-      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-        res |= (*bi).expand(lo,hi);
-      }
-    }
-  }
-  return res;
-}
-
-// cost: O(n^2) in general, but only O(n) for shifting
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi,
-                                   const vect<T,D>& denom) const
-{
-  assert (all(denom > vect<T,D>(0)));
-  bboxset res;
-  {
-    SKIP_NORMALIZE(res);
-    if (all (lo == -hi)) {
-      // Special case for shifting, since this is faster
-      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-        res += (*bi).expand(lo,hi,denom);
-      }
-    } else {
-      // We don't know (yet?) how to shrink a set
-      assert (all ((lo>=0 and hi>=0) or (lo == hi)));
-      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-        res |= (*bi).expand(lo,hi,denom);
-      }
-    }
-  }
-  return res;
-}
-
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::expanded_for (const box& b) const
-{
-  bboxset res;
-  {
-    SKIP_NORMALIZE(res);
-    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-      res |= (*bi).expanded_for(b);
-    }
-  }
-  return res;
-}
-
-// TODO: this is incorrect
-#if 1
-template<typename T, int D>
-bboxset<T,D> bboxset<T,D>::contracted_for (const box& b) const
-{
-  bboxset res;
-  {
-    SKIP_NORMALIZE(res);
-    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
-      res |= (*bi).contracted_for(b);
-    }
-  }
-  return res;
-}
-#endif
-
-
-
-// Equality
-template<typename T, int D>
-bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const
-{
-  return (*this - s).empty();
-}
-
-template<typename T, int D>
-bool bboxset<T,D>::operator< (const bboxset<T,D>& s) const
-{
-  return (*this - s).empty() && not (s - *this).empty();
-}
-
-template<typename T, int D>
-bool bboxset<T,D>::operator>= (const bboxset<T,D>& s) const
-{
-  return s <= *this;
-}
-
-template<typename T, int D>
-bool bboxset<T,D>::operator> (const bboxset<T,D>& s) const
-{
-  return s < *this;
-}
-
-template<typename T, int D>
-bool bboxset<T,D>::operator== (const bboxset<T,D>& s) const
-{
-  return (*this <= s) && (*this >= s);
-}
-
-template<typename T, int D>
-bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const
-{
-  return not (*this == s);
-}
-
-
-
-// Input
-template<typename T,int D>
-istream& bboxset<T,D>::input (istream& is)
-{
-  T Tdummy;
-  try {
-    skipws (is);
-    consume (is, "bboxset<");
-    consume (is, typestring(Tdummy));
-    consume (is, ",");
-    int D_;
-    is >> D_;
-    if (D_ != D) {
-      cout << "Input error: Wrong bboxset dimension " << D_ << ", expected " << D << endl;
-      throw input_error();
-    }
-    consume (is, ">:{");
-    consume (is, "size=");
-    size_type size_;
-    is >> size_;
-    consume (is, ",");
-    consume (is, "setsize=");
-    int setsize_;
-    is >> setsize_;
-    consume (is, ",");
-    consume (is, "set=");
-    is >> bs;
-    consume (is, "}");
-  } catch (input_error & err) {
-    cout << "Input error while reading a bboxset<" << typestring(Tdummy) << "," << D << ">" << endl;
-    throw err;
-  }
-  normalize();
-  return is;
-}
-
-
-
-// Output
-template<typename T,int D>
-ostream& bboxset<T,D>::output (ostream& os) const
-{
-  T Tdummy;
-  os << "bboxset<" << typestring(Tdummy) << "," << D << ">:{"
-     << "size=" << size() << ","
-     << "setsize=" << setsize() << ","
-     << "set=" << bs
-     << "}";
-  return os;
-}
-
-
-
-template class bboxset<int,dim>;
-template size_t memoryof (const bboxset<int,dim>& s);
-template istream& operator>> (istream& is, bboxset<int,dim>& s);
-template ostream& operator<< (ostream& os, const bboxset<int,dim>& s);
-
-#include "dh.hh"
-#include "region.hh"
-
-template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bbox<int,dim> dh::full_dboxes::* const v);
-template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bboxset<int,dim> dh::full_dboxes::* const v);
-template bboxset<int,dim>::bboxset (const vector<region_t>& vb, const bbox<int,dim> region_t::* const v);
diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh
index a8e2ab0ad..463e19f54 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.hh
@@ -1,380 +1,15 @@
 #ifndef BBOXSET_HH
 #define BBOXSET_HH
 
-#include <cassert>
-#include <iostream>
-#include <list>
-#include <set>
-#include <vector>
-
-#include "bbox.hh"
-#include "defs.hh"
-#include "vect.hh"
-
-using namespace std;
-
-
-
-// Choose the implementation of bboxset by #defining exactly one of
-// these
-#undef BBOXSET_SET              // outdated
-#undef BBOXSET_LIST             // well tested
-#define BBOXSET_VECTOR          // brand new
-
-
-
-// Forward declaration
-template<typename T, int D> class bboxset;
-
-// template<typename T,int D>
-// bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2);
-// template<typename T,int D>
-// bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s);
-
-// template<typename T,int D>
-// bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2);
-// template<typename T,int D>
-// bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s);
-
-// Input
-template<typename T,int D>
-istream& operator>> (istream& is, bboxset<T,D>& s);
-
-// Output
-template<typename T,int D>
-ostream& operator<< (ostream& os, const bboxset<T,D>& s);
-
-
-
-// Bounding box set class
-template<typename T, int D>
-class bboxset {
-  
-  // Cost annotations depend on the number of bset elements n, and
-  // assume that normalization is skipped.
-  
-  struct skip_normalize_t {
-    bboxset<T,D>& s;
-    bool const saved_skip_normalize;
-    skip_normalize_t(bboxset<T,D>& s_)
-      : s(s_),
-        saved_skip_normalize(s.skip_normalize)
-    {
-      s.skip_normalize = true;
-    }
-    ~skip_normalize_t()
-    {
-      s.skip_normalize = saved_skip_normalize;
-      s.normalize();
-    }
-  };
-  
-  // Types
-  typedef bbox<T,D> box;
-#ifdef BBOXSET_SET
-  //S typedef set<box> bset;
-#endif
-#ifdef BBOXSET_LIST
-  typedef list<box> bset;
-#endif
-#ifdef BBOXSET_VECTOR
-  typedef vector<box> bset;
-#endif
-  
-  // Fields
-  bset bs;
-  // Invariant:
-  // All bboxes have the same stride.
-  // No bbox is empty.
-  // The bboxes don't overlap.
-  
-  bool skip_normalize;
-  
-public:
-  
-  // Constructors
-  bboxset ();                   // cost: O(1)
-  bboxset (const box& b);       // cost: O(1)
-  bboxset (const bboxset& s);   // cost: O(n)
-  
-  bboxset (const list<box>& lb);
-  bboxset (const vector<box>& lb);
-  bboxset (const vector<list<box> >& vlb);
-  template<typename U>
-  bboxset (const vector<U>& vb, const bbox<T,D> U::* const v);
-  template<typename U>
-  bboxset (const vector<U>& vb, const bboxset U::* const v);
-  
-  static bboxset poison ();
-  
-  // Invariant
-  bool invariant () const CCTK_MEMBER_ATTRIBUTE_PURE;
-  
-private:
-  // Normalisation
-  void normalize ();
-  
-public:
-  // Accessors
-  bool empty () const { return bs.empty(); } // cost: O(1)
-  // T size () const;
-  typedef typename box::size_type size_type;
-  size_type size () const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  int setsize () const { return bs.size(); } // cost: O(1)
-  
-  // Find out whether the bbox contains the point x
-  bool contains (const vect<T,D>& x)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  
-  // Find out whether this bboxset intersects the bbox b
-  bool intersects (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  
-  // Add (bboxes that don't overlap)
-  bboxset& operator+= (const box& b);     // cost: O(1)
-  bboxset& operator+= (const bboxset& s); // cost: O(n)
-  bboxset& add_transfer (bboxset& s);     // cost: O(1)
-  bboxset operator+ (const box& b)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  bboxset operator+ (const bboxset& s)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  
-  // Union
-  bboxset& operator|= (const box& b);     // cost: O(n)
-  bboxset& operator|= (const bboxset& s); // cost: O(n^2)
-  bboxset operator| (const box& b)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  bboxset operator| (const bboxset& s)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n^2)
-  
-  // Intersection
-  bboxset operator& (const box& b)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  bboxset operator& (const bboxset& s)
-    const CCTK_MEMBER_ATTRIBUTE_PURE;     // cost: O(n)
-  bboxset& operator&= (const box& b);     // cost: O(n)
-  bboxset& operator&= (const bboxset& s); // cost: O(n)
-  
-  // Difference
-  // friend bboxset operator- <T,D>(const box& b1, const box& b2);
-  static bboxset minus (const box& b1, const box& b2)
-    CCTK_ATTRIBUTE_PURE;        // cost: O(1)
-  bboxset operator- (const box& b)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  // cost: O(n)
-  bboxset& operator-= (const box& b)
-  {
-    *this = *this - b;
-    assert (invariant());
-    return *this;
-  }
-  bboxset& operator-= (const bboxset& s);     // cost: O(n^2)
-  bboxset operator- (const bboxset& s)
-    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n^2)
-  // friend bboxset operator- <T,D>(const box& b, const bboxset& s);
-  static bboxset minus (const box& b, const bboxset& s)
-    CCTK_ATTRIBUTE_PURE;        // cost: O(n^2)
-  
-  /** Find a bbox containing the whole set.  */
-  box container () const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
-  /** Find the pseudo-inverse.  */
-  bboxset pseudo_inverse (const int n) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  
-  /** Expand (enlarge) the bboxset by multiples of the stride.  */
-  bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi)
-    const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bboxset expand (const vect<vect<T,D>,2>& lohi) const
-  { return expand (lohi[0], lohi[1]); }
-  
-  /** Shift the bboxset by multiples of the stride.  */
-  bboxset shift (const vect<T,D>& v) const
-  { return expand (-v, v); }
-  
-  /** Expand (enlarge) the bboxset by multiples of a fraction of the
-      stride.  */
-  // cost: O(n^2) in general, but only O(n) for shifting
-  bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi,
-                  const vect<T,D>& denom) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  // cost: O(n^2) in general, but only O(n) for shifting
-  bboxset expand (const vect<vect<T,D>,2>& lohi, const vect<T,D>& denom) const
-  { return expand (lohi[0], lohi[1], denom); }
-  
-  /** Shift the bboxset by multiples of a fraction of the stride.  */
-// cost: O(n)
-  bboxset shift (const vect<T,D>& v, const vect<T,D>& denom) const
-  { return expand (-v, v, denom); }
-  
-  /** Find the smallest b-compatible box around this bbox.
-      ("compatible" means having the same stride.)  */
-  bboxset expanded_for (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  
-  // TODO: this is incorrect
-#if 1
-  /** Find the largest b-compatible box inside this bbox.  */
-  bboxset contracted_for (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE;
+#include "bboxset1.hh"
+#ifndef CARPET_NO_BBOXSET2
+#  include "bboxset2.hh"
 #endif
-  
-  // Equality
-  bool operator== (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bool operator!= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bool operator< (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bool operator<= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bool operator> (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  bool operator>= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
-  
-  // Iterators
-  typedef typename bset::const_iterator const_iterator;
-  typedef typename bset::iterator       iterator;
-  
-  const_iterator begin () const { return bs.begin(); }
-  const_iterator end () const   { return bs.end(); }
-//   iterator begin () const { return bs.begin(); }
-//   iterator end () const   { return bs.end(); }
-  
-  // Memory usage
-  size_t memory () const { return memoryof (bs); }
-  
-  // Input
-  istream& input (istream& is);
-  
-  // Output
-  ostream& output (ostream& os) const;
-};
-
-
-
-template<typename T,int D>
-inline bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2) {
-  return bboxset<T,D>(b1) + bboxset<T,D>(b2);
-}
-
-template<typename T,int D>
-inline bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s) {
-  return bboxset<T,D>(b) + s;
-}
-
-// cost: O(1)
-template<typename T,int D>
-inline bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) {
-  return bboxset<T,D>::minus(b1,b2);
-}
-
-// cost: O(n^2)
-template<typename T,int D>
-inline bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) {
-  return bboxset<T,D>::minus(b,s);
-}
-
-
-
-template<typename T,int D>
-inline bboxset<T,D> operator| (const bbox<T,D>& b, const bboxset<T,D>& s) {
-  return s | b;
-}
-
-template<typename T,int D>
-inline bboxset<T,D> operator& (const bbox<T,D>& b, const bboxset<T,D>& s) {
-  return s & b;
-}
-
-
-
-template<typename T,int D>
-inline bool operator== (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) == s;
-}
-
-template<typename T,int D>
-inline bool operator!= (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) != s;
-}
-
-template<typename T,int D>
-inline bool operator< (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) < s;
-}
-
-template<typename T,int D>
-inline bool operator<= (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) <= s;
-}
-
-template<typename T,int D>
-inline bool operator> (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) > s;
-}
-
-template<typename T,int D>
-inline bool operator>= (const bbox<T,D>& b, const bboxset<T,D>& s)
-{
-  return bboxset<T,D>(b) >= s;
-}
-
-
-
-template<typename T,int D>
-inline bool operator== (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s == bboxset<T,D>(b);
-}
-
-template<typename T,int D>
-inline bool operator!= (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s != bboxset<T,D>(b);
-}
-
-template<typename T,int D>
-inline bool operator< (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s < bboxset<T,D>(b);
-}
-
-template<typename T,int D>
-inline bool operator<= (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s <= bboxset<T,D>(b);
-}
-
-template<typename T,int D>
-inline bool operator> (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s > bboxset<T,D>(b);
-}
-
-template<typename T,int D>
-inline bool operator>= (const bboxset<T,D>& s, const bbox<T,D>& b)
-{
-  return s >= bboxset<T,D>(b);
-}
-
-
-
-// Memory usage
-template<typename T, int D>
-inline size_t memoryof (bboxset<T,D> const & s)
-{ return s.memory(); }
-
-
-
-// Input
-template<typename T,int D>
-inline istream& operator>> (istream& is, bboxset<T,D>& s) {
-  return s.input(is);
-}
-
-
-
-// Output
-template<typename T,int D>
-inline ostream& operator<< (ostream& os, const bboxset<T,D>& s) {
-  return s.output(os);
-}
-
 
+#ifdef CARPET_BBOXSET2
+using namespace bboxset2;
+#else
+using namespace bboxset1;
+#endif
 
-#endif // BBOXSET_HH
+#endif  // #ifndef BBOXSET_HH
diff --git a/Carpet/CarpetLib/src/bboxset1.cc b/Carpet/CarpetLib/src/bboxset1.cc
new file mode 100644
index 000000000..def1a6ba6
--- /dev/null
+++ b/Carpet/CarpetLib/src/bboxset1.cc
@@ -0,0 +1,55 @@
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <limits>
+#include <set>
+#include <stack>
+#include <vector>
+
+#include "defs.hh"
+
+#include "bboxset1.hh"
+
+using namespace std;
+
+
+
+namespace bboxset1 {
+
+
+
+template class bboxset<int,1>;
+template void bboxset<int,1>::serialise (set<bbox<int,1> >& s) const;
+template void bboxset<int,1>::serialise (vector<bbox<int,1> >& v) const;
+template size_t memoryof (const bboxset<int,1>& s);
+template istream& operator>> (istream& is, bboxset<int,1>& s);
+template ostream& operator<< (ostream& os, const bboxset<int,1>& s);
+
+template class bboxset<int,2>;
+template void bboxset<int,2>::serialise (set<bbox<int,2> >& s) const;
+template void bboxset<int,2>::serialise (vector<bbox<int,2> >& v) const;
+template size_t memoryof (const bboxset<int,2>& s);
+template istream& operator>> (istream& is, bboxset<int,2>& s);
+template ostream& operator<< (ostream& os, const bboxset<int,2>& s);
+
+template class bboxset<int,3>;
+template void bboxset<int,3>::serialise (set<bbox<int,3> >& s) const;
+template void bboxset<int,3>::serialise (vector<bbox<int,3> >& f) const;
+template size_t memoryof (const bboxset<int,3>& s);
+template istream& operator>> (istream& is, bboxset<int,3>& s);
+template ostream& operator<< (ostream& os, const bboxset<int,3>& s);
+
+} // namespace bboxset1
+
+
+
+#include "dh.hh"
+#include "region.hh"
+
+namespace bboxset1 {
+
+template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bbox<int,dim> dh::full_dboxes::* const v);
+template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bboxset<int,dim> dh::full_dboxes::* const v);
+template bboxset<int,dim>::bboxset (const vector<region_t>& vb, const bbox<int,dim> region_t::* const v);
+
+} // namespace bboxset1
diff --git a/Carpet/CarpetLib/src/bboxset1.hh b/Carpet/CarpetLib/src/bboxset1.hh
new file mode 100644
index 000000000..db760e5ae
--- /dev/null
+++ b/Carpet/CarpetLib/src/bboxset1.hh
@@ -0,0 +1,1189 @@
+#ifndef BBOXSET1_HH
+#define BBOXSET1_HH
+
+#include <cassert>
+#include <iostream>
+#include <list>
+#include <set>
+#include <vector>
+
+#include "bbox.hh"
+#include "defs.hh"
+#include "vect.hh"
+
+using namespace std;
+
+
+
+// Choose the implementation of bboxset by #defining exactly one of
+// these
+#undef BBOXSET_SET              // outdated
+#undef BBOXSET_LIST             // well tested
+#define BBOXSET_VECTOR          // brand new
+
+
+namespace bboxset1 {
+
+
+
+// Forward declaration
+template<typename T, int D> class bboxset;
+
+// template<typename T,int D>
+// bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2);
+// template<typename T,int D>
+// bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s);
+
+// template<typename T,int D>
+// bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2);
+// template<typename T,int D>
+// bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s);
+
+// Input
+template<typename T,int D>
+istream& operator>> (istream& is, bboxset<T,D>& s);
+
+// Output
+template<typename T,int D>
+ostream& operator<< (ostream& os, const bboxset<T,D>& s);
+
+
+
+// Bounding box set class
+template<typename T, int D>
+class bboxset {
+  
+  // Cost annotations depend on the number of bset elements n, and
+  // assume that normalization is skipped.
+  
+  struct skip_normalize_t {
+    bboxset<T,D>& s;
+    bool const saved_skip_normalize;
+    skip_normalize_t(bboxset<T,D>& s_)
+      : s(s_),
+        saved_skip_normalize(s.skip_normalize)
+    {
+      s.skip_normalize = true;
+    }
+    ~skip_normalize_t()
+    {
+      s.skip_normalize = saved_skip_normalize;
+      s.normalize();
+    }
+  };
+  
+  // Types
+  typedef bbox<T,D> box;
+#ifdef BBOXSET_SET
+  //S typedef set<box> bset;
+#endif
+#ifdef BBOXSET_LIST
+  typedef list<box> bset;
+#endif
+#ifdef BBOXSET_VECTOR
+  typedef vector<box> bset;
+#endif
+  
+  // Fields
+  bset bs;
+  // Invariant:
+  // All bboxes have the same stride.
+  // No bbox is empty.
+  // The bboxes don't overlap.
+  
+  bool skip_normalize;
+  
+public:
+  
+  // Constructors
+  bboxset ();                   // cost: O(1)
+  bboxset (const box& b);       // cost: O(1)
+  bboxset (const bboxset& s);   // cost: O(n)
+  
+  bboxset (const list<box>& lb);
+  bboxset (const set<box>& sb);
+  bboxset (const vector<box>& vb);
+  bboxset (const vector<list<box> >& vlb);
+  template<typename U>
+  bboxset (const vector<U>& vb, const bbox<T,D> U::* const v);
+  template<typename U>
+  bboxset (const vector<U>& vb, const bboxset U::* const v);
+  
+  static bboxset poison ();
+  
+  // Invariant
+  bool invariant () const CCTK_MEMBER_ATTRIBUTE_PURE;
+  
+private:
+  // Normalisation
+  void normalize ();
+  
+public:
+  // Accessors
+  bool empty () const { return bs.empty(); } // cost: O(1)
+  // T size () const;
+  typedef typename box::size_type size_type;
+  size_type size () const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  int setsize () const { return bs.size(); } // cost: O(1)
+  
+  // Find out whether the bbox contains the point x
+  bool contains (const vect<T,D>& x)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  
+  // Find out whether this bboxset intersects the bbox b
+  bool intersects (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  
+  // Add (bboxes that don't overlap)
+  bboxset& operator+= (const box& b);     // cost: O(1)
+  bboxset& operator+= (const bboxset& s); // cost: O(n)
+  bboxset& add_transfer (bboxset& s);     // cost: O(1)
+  bboxset operator+ (const box& b)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  bboxset operator+ (const bboxset& s)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  
+  // Union
+  bboxset& operator|= (const box& b);     // cost: O(n)
+  bboxset& operator|= (const bboxset& s); // cost: O(n^2)
+  bboxset operator| (const box& b)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  bboxset operator| (const bboxset& s)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n^2)
+  
+  // Intersection
+  bboxset operator& (const box& b)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  bboxset operator& (const bboxset& s)
+    const CCTK_MEMBER_ATTRIBUTE_PURE;     // cost: O(n)
+  bboxset& operator&= (const box& b);     // cost: O(n)
+  bboxset& operator&= (const bboxset& s); // cost: O(n)
+  
+  // Difference
+  // friend bboxset operator- <T,D>(const box& b1, const box& b2);
+  static bboxset minus (const box& b1, const box& b2)
+    CCTK_ATTRIBUTE_PURE;        // cost: O(1)
+  bboxset operator- (const box& b)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  // cost: O(n)
+  bboxset& operator-= (const box& b)
+  {
+    *this = *this - b;
+    assert (invariant());
+    return *this;
+  }
+  bboxset& operator-= (const bboxset& s);     // cost: O(n^2)
+  bboxset operator- (const bboxset& s)
+    const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n^2)
+  // friend bboxset operator- <T,D>(const box& b, const bboxset& s);
+  static bboxset minus (const box& b, const bboxset& s)
+    CCTK_ATTRIBUTE_PURE;        // cost: O(n^2)
+  
+  /** Find a bbox containing the whole set.  */
+  box container () const CCTK_MEMBER_ATTRIBUTE_PURE; // cost: O(n)
+  /** Find the pseudo-inverse.  */
+  bboxset pseudo_inverse (const int n) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  
+  /** Expand (enlarge) the bboxset by multiples of the stride.  */
+  bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi)
+    const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bboxset expand (const vect<vect<T,D>,2>& lohi) const
+  { return expand (lohi[0], lohi[1]); }
+  
+  /** Shift the bboxset by multiples of the stride.  */
+  bboxset shift (const vect<T,D>& v) const
+  { return expand (-v, v); }
+  
+  /** Expand (enlarge) the bboxset by multiples of a fraction of the
+      stride.  */
+  // cost: O(n^2) in general, but only O(n) for shifting
+  bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi,
+                  const vect<T,D>& denom) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  // cost: O(n^2) in general, but only O(n) for shifting
+  bboxset expand (const vect<vect<T,D>,2>& lohi, const vect<T,D>& denom) const
+  { return expand (lohi[0], lohi[1], denom); }
+  
+  /** Shift the bboxset by multiples of a fraction of the stride.  */
+// cost: O(n)
+  bboxset shift (const vect<T,D>& v, const vect<T,D>& denom) const
+  { return expand (-v, v, denom); }
+  
+  /** Find the smallest b-compatible box around this bbox.
+      ("compatible" means having the same stride.)  */
+  bboxset expanded_for (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  
+  // TODO: this is incorrect
+#if 1
+  /** Find the largest b-compatible box inside this bbox.  */
+  bboxset contracted_for (const box& b) const CCTK_MEMBER_ATTRIBUTE_PURE;
+#endif
+  
+  // Equality
+  bool operator== (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bool operator!= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bool operator< (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bool operator<= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bool operator> (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  bool operator>= (const bboxset& s) const CCTK_MEMBER_ATTRIBUTE_PURE;
+  
+  // Iterators
+  typedef typename bset::const_iterator const_iterator;
+  typedef typename bset::iterator       iterator;
+  
+  const_iterator begin () const { return bs.begin(); }
+  const_iterator end () const   { return bs.end(); }
+//   iterator begin () const { return bs.begin(); }
+//   iterator end () const   { return bs.end(); }
+  
+  template<typename C>
+  void serialise (C& out) const;
+  
+  // Memory usage
+  size_t memory () const { return memoryof (bs); }
+  
+  // Input
+  istream& input (istream& is);
+  
+  // Output
+  ostream& output (ostream& os) const;
+};
+
+
+
+template<typename T,int D>
+inline bboxset<T,D> operator+ (const bbox<T,D>& b1, const bbox<T,D>& b2) {
+  return bboxset<T,D>(b1) + bboxset<T,D>(b2);
+}
+
+template<typename T,int D>
+inline bboxset<T,D> operator+ (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return bboxset<T,D>(b) + s;
+}
+
+// cost: O(1)
+template<typename T,int D>
+inline bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) {
+  return bboxset<T,D>::minus(b1,b2);
+}
+
+// cost: O(n^2)
+template<typename T,int D>
+inline bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return bboxset<T,D>::minus(b,s);
+}
+
+
+
+template<typename T,int D>
+inline bboxset<T,D> operator| (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return s | b;
+}
+
+template<typename T,int D>
+inline bboxset<T,D> operator& (const bbox<T,D>& b, const bboxset<T,D>& s) {
+  return s & b;
+}
+
+
+
+template<typename T,int D>
+inline bool operator== (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) == s;
+}
+
+template<typename T,int D>
+inline bool operator!= (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) != s;
+}
+
+template<typename T,int D>
+inline bool operator< (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) < s;
+}
+
+template<typename T,int D>
+inline bool operator<= (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) <= s;
+}
+
+template<typename T,int D>
+inline bool operator> (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) > s;
+}
+
+template<typename T,int D>
+inline bool operator>= (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  return bboxset<T,D>(b) >= s;
+}
+
+
+
+template<typename T,int D>
+inline bool operator== (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s == bboxset<T,D>(b);
+}
+
+template<typename T,int D>
+inline bool operator!= (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s != bboxset<T,D>(b);
+}
+
+template<typename T,int D>
+inline bool operator< (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s < bboxset<T,D>(b);
+}
+
+template<typename T,int D>
+inline bool operator<= (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s <= bboxset<T,D>(b);
+}
+
+template<typename T,int D>
+inline bool operator> (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s > bboxset<T,D>(b);
+}
+
+template<typename T,int D>
+inline bool operator>= (const bboxset<T,D>& s, const bbox<T,D>& b)
+{
+  return s >= bboxset<T,D>(b);
+}
+
+
+
+// Memory usage
+template<typename T, int D>
+inline size_t memoryof (bboxset<T,D> const & s)
+{ return s.memory(); }
+
+
+
+// Input
+template<typename T,int D>
+inline istream& operator>> (istream& is, bboxset<T,D>& s) {
+  return s.input(is);
+}
+
+
+
+// Output
+template<typename T,int D>
+inline ostream& operator<< (ostream& os, const bboxset<T,D>& s) {
+  return s.output(os);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+//  Implementation
+////////////////////////////////////////////////////////////////////////////////
+
+#define SKIP_NORMALIZE(s) skip_normalize_t skip_normalize_v(s)
+// #define SKIP_NORMALIZE(s) do { } while(0)
+
+// Constructors
+template<typename T, int D>
+bboxset<T,D>::bboxset ()
+  : skip_normalize(false)
+{
+  assert (invariant());
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const box& b)
+  : skip_normalize(false)
+{
+  //S if (not b.empty()) bs.insert(b);
+  if (not b.empty()) bs.push_back(b);
+  assert (invariant());
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const bboxset& s)
+  : bs(s.bs), skip_normalize(false)
+{
+  assert (invariant());
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const list<box>& lb)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename list<box>::const_iterator
+         li = lb.begin(), le = lb.end(); li != le; ++ li)
+  {
+    *this |= *li;
+  }
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const set<box>& sb)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename set<box>::const_iterator
+         vi = sb.begin(), ve = sb.end(); vi != ve; ++ vi)
+  {
+    *this |= *vi;
+  }
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const vector<box>& vb)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename vector<box>::const_iterator
+         vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi)
+  {
+    *this |= *vi;
+  }
+}
+
+template<typename T, int D>
+bboxset<T,D>::bboxset (const vector<list<box> >& vlb)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename vector<list<box> >::const_iterator
+         vli = vlb.begin(), vle = vlb.end(); vli != vle; ++ vli)
+  {
+    *this |= bboxset (*vli);
+  }
+}
+
+template<typename T, int D>
+template<typename U>
+bboxset<T,D>::bboxset (const vector<U>& vb, const bbox<T,D> U::* const v)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename vector<U>::const_iterator
+         vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi)
+  {
+    *this |= (*vi).*v;
+  }
+}
+
+template<typename T, int D>
+template<typename U>
+bboxset<T,D>::bboxset (const vector<U>& vb, const bboxset U::* const v)
+  : skip_normalize(false)
+{
+  SKIP_NORMALIZE(*this);
+  for (typename vector<U>::const_iterator
+         vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi)
+  {
+    *this |= (*vi).*v;
+  }
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::poison () {
+  return bboxset (bbox<T,D>::poison());
+}
+
+
+
+// Invariant
+template<typename T, int D>
+bool bboxset<T,D>::invariant () const {
+  // This is very slow when there are many bboxes
+#if 0 && defined(CARPET_DEBUG)
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    if ((*bi).empty()) return false;
+    if (not (*bi).is_aligned_with(*bs.begin())) return false;
+    // check for overlap (quadratic -- expensive)
+    for (const_iterator bi2=begin(); bi2!=bi; ++bi2) {
+      if ((*bi2).intersects(*bi)) return false;
+    }
+  }
+#endif
+  return true;
+}
+
+
+
+// Normalisation
+template<typename T, int D>
+void bboxset<T,D>::normalize ()
+{
+  if (skip_normalize) return;
+  assert (invariant());
+  
+  bboxset const oldbs = * this;
+  size_type const oldsize = this->size();
+  
+  // Split all bboxes into small pieces which have all their
+  // boundaries aligned.
+  for (int d=0; d<D; ++d) {
+    // Find all boundaries
+    //S typedef set<T> buf;
+    typedef vector<T> buf;
+    buf sbnds;
+    sbnds.reserve (2 * bs.size());
+    for (typename bset::const_iterator
+           si = bs.begin(), se = bs.end(); si != se; ++ si)
+    {
+      box const & b = * si;
+      int const bstr = b.stride()[d];
+      int const blo = b.lower()[d];
+      int const bhi = b.upper()[d] + bstr;
+      //S sbnds.insert (blo);
+      //S sbnds.insert (bhi);
+      sbnds.push_back (blo);
+      sbnds.push_back (bhi);
+    }
+    sort (sbnds.begin(), sbnds.end());
+    typename buf::iterator const last = unique (sbnds.begin(), sbnds.end());
+    sbnds.resize (last - sbnds.begin());
+    // Split bboxes
+    bset nbs;
+    for (typename bset::const_iterator
+           si = bs.begin(), se = bs.end(); si != se; ++ si)
+    {
+      box const & b = * si;
+      int const bstr = b.stride()[d];
+      int const blo = b.lower()[d];
+      int const bhi = b.upper()[d] + bstr;
+      typename buf::const_iterator const ilo
+        = lower_bound (sbnds.begin(), sbnds.end(), blo);
+      typename buf::const_iterator const ihi
+        = lower_bound (sbnds.begin(), sbnds.end(), bhi);
+      assert (ilo != sbnds.end());
+      assert (ihi != sbnds.end());
+      assert (* ilo == blo);
+      assert (* ihi == bhi);
+      // Split one bbox
+      for (typename buf::const_iterator curr = ilo; curr != ihi; ++ curr) {
+        typename buf::const_iterator next = curr;
+        advance (next, 1);
+        int const nblo = * curr;
+        int const nbhi = * next;
+        assert (nbhi > nblo);   // ensure that the set remains sorted
+        box const nb (b.lower().replace(d, nblo),
+                      b.upper().replace(d, nbhi - bstr),
+                      b.stride());
+        //S nbs.insert (nb);
+        nbs.push_back (nb);
+      }
+    }
+    // Replace old set
+    bs.clear ();
+    bs.swap (nbs);
+    assert (invariant());
+  }
+  
+  // Combine bboxes if possible
+  for (int d=0; d<D; ++d) {
+    bset nbs;
+    while (not bs.empty()) {
+      typename bset::iterator si = bs.begin();
+      assert (si != bs.end());
+      
+      box const b = * si;
+      int const bstr = b.stride()[d];
+      int const blo = b.lower()[d];
+      int const bhi = b.upper()[d] + bstr;
+      
+      for (typename bset::iterator
+             nsi = nbs.begin(), nse = nbs.end(); nsi != nse; ++ nsi)
+      {
+        box const nb = * nsi;
+        int const nblo = nb.lower()[d];
+        int const nbhi = nb.upper()[d] + bstr;
+        
+        box const mb (nb.lower().replace(d, blo),
+                      nb.upper().replace(d, bhi - bstr),
+                      nb.stride());
+        
+        // Check whether the other dimensions match
+        if (b == mb) {
+          // Check whether the bboxes are adjacent in this dimension
+          if (nbhi == blo) {
+            // Combine boxes, nb < b
+            box const cb (b.lower().replace(d, nblo),
+                          b.upper(),
+                          b.stride());
+            bs.erase (si);
+            nbs.erase (nsi);
+            //S bs.insert (cb);
+            bs.push_back (cb);
+            goto done;
+          } else if (bhi == nblo) {
+            // Combine boxes, b < nb
+            box const cb (b.lower(),
+                          b.upper().replace(d, nbhi - bstr),
+                          b.stride());
+            bs.erase (si);
+            nbs.erase (nsi);
+            //S bs.insert (cb);
+            bs.push_back (cb);
+            goto done;
+          }
+        }
+      }
+      bs.erase (si);
+      //S nbs.insert (b);
+      nbs.push_back (b);
+    done:;
+    }
+    bs.swap (nbs);
+    assert (invariant());
+  }
+  
+  size_type const newsize = this->size();
+  
+  // Can't use operators on *this since these would call normalize again
+  // assert (*this == oldbs);
+  assert (newsize == oldsize);
+}
+
+
+
+// Accessors
+// cost: O(n)
+template<typename T, int D>
+typename bboxset<T,D>::size_type bboxset<T,D>::size () const
+{
+  size_type s = 0;
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    const size_type bsz = (*bi).size();
+    assert (numeric_limits<size_type>::max() - bsz >= s);
+    s += bsz;
+  }
+  return s;
+}
+
+
+
+// Queries
+
+// Containment
+// cost: O(n)
+template<typename T, int D>
+bool bboxset<T,D>::contains (const vect<T,D>& x) const
+{
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    if ((*bi).contains(x)) return true;
+  }
+  return false;
+}
+
+// Intersection
+// cost: O(n)
+template<typename T, int D>
+bool bboxset<T,D>::intersects (const box& b) const
+{
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    if ((*bi).intersects(b)) return true;
+  }
+  return false;
+}
+
+
+
+// Add (bboxes that don't overlap)
+// cost: O(1)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator+= (const box& b)
+{
+  if (b.empty()) return *this;
+  // This is very slow when there are many bboxes
+#if 0 && defined(CARPET_DEBUG)
+  // check for overlap
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    assert (not (*bi).intersects(b));
+  }
+#endif
+  //S bs.insert(b);
+  bs.push_back(b);
+  assert (invariant());
+  normalize();
+  return *this;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s)
+{
+  SKIP_NORMALIZE(*this);
+  for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
+    *this += *bi;
+  }
+  assert (invariant());
+  return *this;
+}
+
+//L cost: O(1)
+// cost O(n)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::add_transfer (bboxset& s)
+{
+#ifdef BBOXSET_LIST
+  bs.splice (bs.end(), s.bs);
+#else
+  bs.insert (bs.end(), s.bs.begin(), s.bs.end());
+  s.bs.clear ();
+#endif
+  assert (invariant());
+  normalize();
+  return *this;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator+ (const box& b) const
+{
+  bboxset r(*this);
+  r += b;
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator+ (const bboxset& s) const
+{
+  bboxset r(*this);
+  r += s;
+  assert (r.invariant());
+  return r;
+}
+
+
+
+// Union
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator|= (const box& b)
+{
+  if (b.empty()) return *this;
+#if 0
+  // this has a cost of O(n^2)
+  bboxset tmp = b - *this;
+  add_transfer (tmp);
+#else
+  // this has a cost of O(n)
+  bset oldbs;
+  oldbs.swap(bs);
+  bs.push_back(b);
+  SKIP_NORMALIZE(*this);
+  for (const_iterator bi=oldbs.begin(), be=oldbs.end(); bi!=be; ++bi) {
+    bboxset tmp = *bi - b;
+    add_transfer(tmp);
+  }
+#endif
+  assert (invariant());
+  return *this;
+}
+
+// cost: O(n^2)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s)
+{
+  bboxset tmp = s - *this;
+  add_transfer (tmp);
+  assert (invariant());
+  return *this;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator| (const box& b) const
+{
+  bboxset r(*this);
+  r |= b;
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n^2)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator| (const bboxset& s) const
+{
+  bboxset r(*this);
+  r |= s;
+  assert (r.invariant());
+  return r;
+}
+
+
+
+// Intersection
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator& (const box& b) const
+{
+  // start with an empty set
+  bboxset r;
+  {
+    SKIP_NORMALIZE(r);
+    // walk all my elements
+    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+      // insert the intersection with the bbox
+      r += *bi & b;
+    }
+  }
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const
+{
+  // start with an empty set
+  bboxset r;
+  {
+    SKIP_NORMALIZE(r);
+    // walk all the bboxes
+    for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
+      // insert the intersection with this bbox
+      bboxset tmp = *this & *bi;
+      r.add_transfer (tmp);
+    }
+  }
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator&= (const box& b)
+{
+  *this = *this & b;
+  assert (invariant());
+  return *this;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s)
+{
+  *this = *this & s;
+  assert (invariant());
+  return *this;
+}
+
+
+
+// Difference
+// cost: O(1)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+  assert (b1.is_aligned_with(b2));
+  if (b1.empty()) return bboxset<T,D>();
+  if (not b1.intersects(b2)) return bboxset<T,D>(b1);
+  vect<T,D> const b2lo = max (b1.lower(), b2.lower());
+  vect<T,D> const b2up = min (b1.upper(), b2.upper());
+  vect<T,D> const & b1lo = b1.lower();
+  vect<T,D> const & b1up = b1.upper();
+  vect<T,D> const & str = b1.stride();
+  bboxset<T,D> r;
+  SKIP_NORMALIZE(r);
+  {
+    for (int d=0; d<D; ++d) {
+      // make resulting bboxes as large as possible in x-direction
+      // (for better consumption by Fortranly ordered arrays)
+      vect<T,D> lb, ub;
+      bbox<T,D> b;
+      for (int dd=0; dd<D; ++dd) {
+        if (dd<d) {
+          lb[dd] = b2lo[dd];
+          ub[dd] = b2up[dd];
+        } else if (dd>d) {
+          lb[dd] = b1lo[dd];
+          ub[dd] = b1up[dd];
+        }
+      }
+      lb[d] = b1lo[d];
+      ub[d] = b2lo[d] - str[d];
+      b = bbox<T,D>(lb,ub,str);
+      r += b;
+      lb[d] = b2up[d] + str[d];
+      ub[d] = b1up[d];
+      b = bbox<T,D>(lb,ub,str);
+      r += b;
+    }
+  }
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator- (const box& b) const
+{
+  // start with an empty set
+  bboxset r;
+  {
+    SKIP_NORMALIZE(r);
+    // walk all my elements
+    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+      // insert the difference with the bbox
+      bboxset tmp = *bi - b;
+      r.add_transfer (tmp);
+    }
+  }
+  assert (r.invariant());
+  return r;
+}
+  
+// cost: O(n^2)
+template<typename T, int D>
+bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s)
+{
+  SKIP_NORMALIZE(*this);
+  for (const_iterator bi=s.begin(), be=s.end(); bi!=be; ++bi) {
+    *this -= *bi;
+  }
+  assert (invariant());
+  return *this;
+}
+
+// cost: O(n^2)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const
+{
+  bboxset r(*this);
+  r -= s;
+  assert (r.invariant());
+  return r;
+}
+
+// cost: O(n^2)
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s)
+{
+  bboxset<T,D> r = bboxset<T,D>(b) - s;
+  assert (r.invariant());
+  return r;
+}
+
+
+
+// cost: O(n)
+template<typename T, int D>
+typename bboxset<T,D>::box bboxset<T,D>::container () const
+{
+  box b;
+  for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+    b = b.expanded_containing(*bi);
+  }
+  return b;
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::pseudo_inverse (const int n) const
+{
+  assert (not empty());
+  box const c = container().expand(n,n);
+  return c - *this;
+}
+
+// cost: O(n^2) in general, but only O(n) for shifting
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi)
+  const
+{
+  bboxset res;
+  {
+    SKIP_NORMALIZE(res);
+    if (all (lo == -hi)) {
+      // Special case for shifting, since this is faster
+      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+        res += (*bi).expand(lo,hi);
+      }
+    } else {
+      // We don't know (yet?) how to shrink a set
+      assert (all (lo>=0 and hi>=0));
+      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+        res |= (*bi).expand(lo,hi);
+      }
+    }
+  }
+  return res;
+}
+
+// cost: O(n^2) in general, but only O(n) for shifting
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi,
+                                   const vect<T,D>& denom) const
+{
+  assert (all(denom > vect<T,D>(0)));
+  bboxset res;
+  {
+    SKIP_NORMALIZE(res);
+    if (all (lo == -hi)) {
+      // Special case for shifting, since this is faster
+      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+        res += (*bi).expand(lo,hi,denom);
+      }
+    } else {
+      // We don't know (yet?) how to shrink a set
+      assert (all ((lo>=0 and hi>=0) or (lo == hi)));
+      for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+        res |= (*bi).expand(lo,hi,denom);
+      }
+    }
+  }
+  return res;
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::expanded_for (const box& b) const
+{
+  bboxset res;
+  {
+    SKIP_NORMALIZE(res);
+    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+      res |= (*bi).expanded_for(b);
+    }
+  }
+  return res;
+}
+
+// TODO: this is incorrect
+#if 1
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::contracted_for (const box& b) const
+{
+  bboxset res;
+  {
+    SKIP_NORMALIZE(res);
+    for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) {
+      res |= (*bi).contracted_for(b);
+    }
+  }
+  return res;
+}
+#endif
+
+
+
+// Equality
+template<typename T, int D>
+bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const
+{
+  return (*this - s).empty();
+}
+
+template<typename T, int D>
+bool bboxset<T,D>::operator< (const bboxset<T,D>& s) const
+{
+  return (*this - s).empty() && not (s - *this).empty();
+}
+
+template<typename T, int D>
+bool bboxset<T,D>::operator>= (const bboxset<T,D>& s) const
+{
+  return s <= *this;
+}
+
+template<typename T, int D>
+bool bboxset<T,D>::operator> (const bboxset<T,D>& s) const
+{
+  return s < *this;
+}
+
+template<typename T, int D>
+bool bboxset<T,D>::operator== (const bboxset<T,D>& s) const
+{
+  return (*this <= s) && (*this >= s);
+}
+
+template<typename T, int D>
+bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const
+{
+  return not (*this == s);
+}
+
+
+
+// Serialise
+template<typename T, int D>
+template<typename C>
+void bboxset<T,D>::serialise (C& out) const
+{
+  for (bboxset::const_iterator it = begin(); it != end(); ++it) {
+    out.insert(out.end(), *it);
+  }
+}
+
+
+
+// Input
+template<typename T,int D>
+istream& bboxset<T,D>::input (istream& is)
+{
+  T Tdummy;
+  try {
+    skipws (is);
+    consume (is, "bboxset<");
+    consume (is, typestring(Tdummy));
+    consume (is, ",");
+    int D_;
+    is >> D_;
+    if (D_ != D) {
+      cout << "Input error: Wrong bboxset dimension " << D_ << ", expected " << D << endl;
+      throw input_error();
+    }
+    consume (is, ">:{");
+    consume (is, "size=");
+    size_type size_;
+    is >> size_;
+    consume (is, ",");
+    consume (is, "setsize=");
+    int setsize_;
+    is >> setsize_;
+    consume (is, ",");
+    consume (is, "set=");
+    is >> bs;
+    consume (is, "}");
+  } catch (input_error & err) {
+    cout << "Input error while reading a bboxset<" << typestring(Tdummy) << "," << D << ">" << endl;
+    throw err;
+  }
+  normalize();
+  return is;
+}
+
+
+
+// Output
+template<typename T,int D>
+ostream& bboxset<T,D>::output (ostream& os) const
+{
+  T Tdummy;
+  os << "bboxset<" << typestring(Tdummy) << "," << D << ">:{"
+     << "size=" << size() << ","
+     << "setsize=" << setsize() << ","
+     << "set=" << bs
+     << "}";
+  return os;
+}
+
+
+
+#undef SKIP_NORMALIZE
+
+
+
+} // namespace bboxset1
+
+
+
+#endif // BBOXSET1_HH
diff --git a/Carpet/CarpetLib/src/bboxset2.cc b/Carpet/CarpetLib/src/bboxset2.cc
new file mode 100644
index 000000000..447d04951
--- /dev/null
+++ b/Carpet/CarpetLib/src/bboxset2.cc
@@ -0,0 +1,38 @@
+#include <cctk.h>
+
+#ifndef CARPET_NO_BBOXSET2
+
+#include "bboxset2.hh"
+
+
+
+namespace bboxset2 {
+  
+  template class bboxset<int,0>;
+  
+  template class bboxset<int,1>;
+  template void bboxset<int,1>::serialise(set<bbox>&) const;
+  
+  template class bboxset<int,2>;
+  template void bboxset<int,2>::serialise(set<bbox>&) const;
+  
+  template class bboxset<int,3>;
+  template void bboxset<int,3>::serialise(set<bbox>&) const;
+  template void bboxset<int,3>::serialise(vector<bbox>&) const;
+  
+} // namespace bboxset2
+
+
+
+#include "dh.hh"
+#include "region.hh"
+
+namespace bboxset2 {
+  
+  template bboxset<int,3>::bboxset(const vector<dh::full_dboxes>&, const bbox dh::full_dboxes::*);
+  template bboxset<int,3>::bboxset(const vector<dh::full_dboxes>&, const bboxset dh::full_dboxes::*);
+  template bboxset<int,3>::bboxset(const vector<region_t>&, const bbox region_t::*);
+  
+} // namespace bboxset2
+
+#endif  // #ifndef CARPET_NO_BBOXSET2
diff --git a/Carpet/CarpetLib/src/bboxset2.hh b/Carpet/CarpetLib/src/bboxset2.hh
new file mode 100644
index 000000000..bac295377
--- /dev/null
+++ b/Carpet/CarpetLib/src/bboxset2.hh
@@ -0,0 +1,1369 @@
+#ifndef BBOXSET2_HH
+#define BBOXSET2_HH
+
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+#include <map>
+#include <memory>
+#include <set>
+#include <sstream>
+#include <utility>
+#include <vector>
+
+#include <cctk.h>
+
+#include "bbox.hh"
+#include "defs.hh"
+#include "vect.hh"
+
+using namespace std;
+
+
+
+#ifndef CARPET_NO_BBOXSET2
+
+namespace bboxset2 {
+
+
+
+template<typename T, int D>
+class bboxset {
+  template<typename, int> friend class bboxset;
+  typedef ::vect<T,D> vect;
+  typedef ::bbox<T,D> bbox;
+  typedef ::vect<T,D-1> vect1;
+  typedef ::bbox<T,D-1> bbox1;
+  typedef bboxset<T,D-1> bboxset1;
+  
+#if 0 && !defined __ICC
+  template<typename U>
+  using ptr = unique_ptr<U>;
+#else
+  template<typename U>
+  using ptr = shared_ptr<U>;
+#endif
+  
+  typedef map<T, ptr<bboxset1>> subsets_t;
+  subsets_t subsets;
+  
+  vect stride, offset;
+  
+  bool is_poison_;
+  
+  template<typename F>
+  void traverse_subsets(const F& f) const
+  {
+    bboxset1 decoded_subset;
+    assert(decoded_subset.empty());
+    for (const auto& pos_subset: subsets) {
+      const T& pos = pos_subset.first;
+      const bboxset1* const subsetp = pos_subset.second.get();
+      decoded_subset ^= *subsetp;
+      f(pos, decoded_subset);
+    }
+    assert(decoded_subset.empty());
+  }
+  
+  template<typename F>
+  void traverse_subsets(const F& f, const bboxset& other) const
+  {
+    bboxset1 decoded_subset0;
+    bboxset1 decoded_subset1;
+    assert(decoded_subset0.empty());
+    assert(decoded_subset1.empty());
+    
+    typedef typename subsets_t::const_iterator subsets_iter_t;
+    subsets_iter_t iter0 = subsets.begin();
+    subsets_iter_t iter1 = other.subsets.begin();
+    subsets_iter_t const end0 = subsets.end();
+    subsets_iter_t const end1 = other.subsets.end();
+    while (iter0!=end0 or iter1!=end1) {
+      const T next_pos0 = iter0!=end0 ? iter0->first : numeric_limits<T>::max();
+      const T next_pos1 = iter1!=end1 ? iter1->first : numeric_limits<T>::max();
+      const T pos = min(next_pos0, next_pos1);
+      const bool active0 = next_pos0==pos;
+      const bool active1 = next_pos1==pos;
+      const bboxset1* const subset0p = active0 ? iter0->second.get() : 0;
+      const bboxset1* const subset1p = active1 ? iter1->second.get() : 0;
+      if (active0) decoded_subset0 ^= *subset0p;
+      if (active1) decoded_subset1 ^= *subset1p;
+      
+      f(pos, decoded_subset0, decoded_subset1);
+      
+      if (active0) ++iter0;
+      if (active1) ++iter1;
+    }
+    assert(decoded_subset0.empty());
+    assert(decoded_subset1.empty());
+  }
+  
+  template<typename F>
+  bboxset binary_operator(const F& op, const bboxset& other) const
+  {
+    assert(not is_poison() and not other.is_poison());
+    
+#if 0
+    // TODO: This assumes that the empty set is a neutral element for
+    // the operator
+    if (other.empty()) return *this;
+    if (empty()) return other;
+    assert(all(other.stride == stride));
+    assert(all(offset == other.offset));
+#endif
+    
+    if (not empty() and not other.empty()) {
+      assert(all(stride == other.stride));
+      // TODO
+      if (not (all(offset == other.offset))) {
+        cout << "*this=" << *this << "\n"
+             << "other=" << other << "\n";
+      }
+      assert(all(offset == other.offset));
+    }
+    
+    bboxset res;
+    if (empty()) {
+      res.stride = other.stride;
+      res.offset = other.offset;
+    } else {
+      res.stride = stride;
+      res.offset = offset;
+    }
+    
+    bboxset1 old_decoded_subsetr;
+    traverse_subsets
+      ([&](const T& pos,
+           const bboxset1& decoded_subset0,
+           const bboxset1& decoded_subset1)
+       {
+         bboxset1 decoded_subsetr = op(decoded_subset0, decoded_subset1);
+         const auto subsetrp =
+           make_shared<bboxset1>(decoded_subsetr ^ old_decoded_subsetr);
+         if (not subsetrp->empty()) {
+           res.subsets.insert(res.subsets.end(), make_pair(pos, subsetrp));
+         }
+         swap(old_decoded_subsetr, decoded_subsetr);
+       },
+       other);
+    assert(old_decoded_subsetr.empty());
+    
+    return res;
+  }
+  
+public:
+  
+  /** Invariant */
+  bool invariant() const
+  {
+    if (is_poison() and not empty()) return false;
+    if (any(stride <= vect(0))) return false;
+    if (any(offset < vect(0) or offset >= stride)) return false;
+    for (const auto& pos_subset: subsets) {
+      if (pos_subset.second.get()->empty()) return false;
+      if (any(pos_subset.second.get()->stride != init(stride))) return false;
+      if (any(pos_subset.second.get()->offset != init(offset))) return false;
+    }
+    if (chi_size() % 2 != 0) return false;
+    return true;
+  }
+  
+  
+  
+  /** Create empty set */
+  bboxset(): stride(vect(1)), offset(vect(0)), is_poison_(false)
+  {
+  }
+  
+  /** Return a is_poisoned set */
+  static bboxset poison()
+  {
+    bboxset bs;
+    bs.is_poison_ = true;
+    return bs;
+  }
+  
+  /** Copy constructor */
+  bboxset(const bboxset& other);
+  
+  /** Copy constructor */
+  // bboxset(bboxset&& other);
+  
+  /** Assignment */
+  bboxset& operator=(const bboxset& other);
+  
+  /** Assignment */
+  // bboxset& operator=(bboxset&& other);
+  
+  /** Create set from bbox */
+  bboxset(const bbox& b);
+  
+  /** Create set from container of bboxes, bboxsets, or container
+      thereof */
+  template<typename C>
+  bboxset(const C& elts);
+  
+  /** Create set from container of structs containing a bbox, bboxset,
+      or container thereof */
+  template<typename C, typename S, typename B>
+  bboxset(const C& elts, const B S::* const mptr);
+  
+  
+  
+  /** Test for is_poison */
+  bool is_poison() const
+  {
+    return is_poison_;
+  }
+  
+  /** Test for emptiness */
+  bool empty() const
+  {
+    assert(not is_poison());
+    return subsets.empty();
+  }
+  
+  /** Find a point contained in this set */
+  vect front() const;
+  
+  /** Number of characteristic points */
+  T chi_size() const;
+  
+  /** Number of elements */
+  T size() const;
+  
+  /** Container (min and max) */
+  bbox container() const;
+  
+  
+  
+  /** Test for equality */
+  bool operator==(const bboxset& other) const;
+  
+  /** Test for is-subset-of */
+  bool operator<=(const bboxset& other) const;
+  
+  /** Test for is-strict-subset-of */
+  bool operator<(const bboxset& other) const;
+  
+  bool operator!=(const bboxset& other) const
+  {
+    return not (*this == other);
+  }
+  
+  bool operator>=(const bboxset& other) const
+  {
+    return other <= *this;
+  }
+  
+  bool operator>(const bboxset& other) const
+  {
+    return other < *this;
+  }
+  
+  bool contains(const vect& v) const
+  {
+    return bbox(v, v, stride) <= *this;
+  }
+  
+  bool intersects(const bbox& b) const
+  {
+    return not (b & *this).empty();
+  }
+  
+  
+  
+  /** Symmetric set difference */
+  bboxset operator^(const bboxset& other) const;
+  
+  /** Symmetric set difference */
+  bboxset& operator^=(const bboxset& other);
+  
+  /** Set intersection */
+  bboxset operator&(const bboxset& other) const;
+  
+  /** Set intersection */
+  bboxset& operator&=(const bboxset& other);
+  
+  /** Set Union */
+  bboxset operator|(const bboxset& other) const;
+  
+  /** Set union */
+  bboxset& operator|=(const bboxset& other);
+  
+  /** Symmetric set union */
+  bboxset operator+(const bboxset& other) const;
+  
+  /** Symmetric set union */
+  bboxset& operator+=(const bboxset& other);
+  
+  /** Set difference */
+  bboxset operator-(const bboxset& other) const;
+  
+  /** Set difference */
+  bboxset& operator-=(const bboxset& other);
+  
+  
+  
+  /** Operators taking bboxes as arguments */
+  bboxset& operator&=(const bbox& b)
+  {
+    return *this &= bboxset(b);
+  }
+  bboxset& operator^=(const bbox& b)
+  {
+    return *this ^= bboxset(b);
+  }
+  bboxset& operator|=(const bbox& b)
+  {
+    return *this |= bboxset(b);
+  }
+  bboxset& operator+=(const bbox& b)
+  {
+    return *this += bboxset(b);
+  }
+  bboxset& operator-=(const bbox& b)
+  {
+    return *this -= bboxset(b);
+  }
+  
+  
+  
+  /** Shift all points */
+  bboxset shift(const vect& dist, const vect& dist_denom = vect(1)) const;
+  
+  /** Expand the set (convolute with a bbox) */
+  bboxset expand(const vect& lo, const vect& hi) const;
+  
+  /** Expand the set (convolute with a bbox) */
+  bboxset expand(const ::vect<vect,2>& lohi) const
+  {
+    return expand(lohi[0], lohi[1]);
+  }
+  
+  
+  
+private:
+  static bool strides_are_compatible(const vect& str1, const vect& str2);
+public:
+  
+  /** Expande the set (changing the stride): find the smallest
+      b-compatible bboxset around this bboxset */
+  // Possible implementation:
+  // - convert bboxset into set of bboxes
+  // - expand each bbox
+  // - create new bboxset as union of these bboxes
+  // Alternative, possibly more efficient implementation:
+  // - for each dimension:
+  //   - BS1 := shifted left to align with bbox
+  //   - BS2 := shifted right to align with bbox
+  //   - create union of BS1 and BS2
+  bboxset expanded_for(const bbox& target) const;
+  
+  /** Shrink the set (changing the stride): find the largest
+      b-compatible bboxset inside this bboxset */
+  bboxset contracted_for(const bbox& target) const;
+  
+  /** Expande the set (changing the stride): find the smallest open
+      b-compatible bboxset around this bboxset */
+  bboxset anti_contracted_for(const bbox& target) const;
+  
+  
+  
+  /** Serialise the set */
+  template<typename C>
+  void serialise(C& out) const;
+  
+  /** Iterate over a serialised set */
+private:
+  typedef vector<bbox> iter_memo_t;
+public:
+  typedef typename iter_memo_t::const_iterator const_iterator;
+private:
+  mutable iter_memo_t iter_memo;
+public:
+  int setsize() const
+  {
+    iter_memo_t im;
+    serialise(im);
+    return im.size();
+  }
+  const_iterator begin() const
+  {
+    iter_memo.clear();
+    serialise(iter_memo);
+    return iter_memo.begin();
+  }
+  const_iterator end() const
+  {
+    return iter_memo.end();
+  }
+  
+  
+  /** Memory usage */
+  size_t memory() const;
+  
+  /** Input */
+  istream& input(istream& is);
+  
+  /** Output */
+  ostream& output(ostream& os) const;
+};
+
+
+
+template<typename T>
+class bboxset<T,0> {
+  template<typename, int> friend class bboxset;
+  typedef ::vect<T,0> vect;
+  typedef ::bbox<T,0> bbox;
+  
+  bool state;
+  vect stride, offset;
+  bool is_poison_;
+public:
+  bool invariant() const
+  {
+    if (is_poison() and state) return false;
+    return true;
+  }
+  
+  bboxset():
+    state(false), is_poison_(false)
+  {
+  }
+  static bboxset poison()
+  {
+    bboxset bs; bs.is_poison_ = true; return bs;
+  }
+  bboxset(const bbox& b):
+    state(true), is_poison_(false)
+  {
+    assert(not b.is_poison());
+  }
+  bboxset(const bboxset& other):
+    state(other.state), is_poison_(other.is_poison_)
+  {
+  }
+  bboxset& operator=(const bboxset& other)
+  {
+    state = other.state;
+    is_poison_ = other.is_poison_;
+    return *this;
+  }
+  
+  bool is_poison() const
+  {
+    return is_poison_;
+  }
+  bool empty() const
+  {
+    assert(not is_poison());
+    return not state;
+  }
+  vect front() const
+  {
+    assert(not is_poison());
+    assert(not empty());
+    return vect();
+  }
+  T chi_size() const
+  {
+    assert(not is_poison());
+    return 1;
+  }
+  T size() const
+  {
+    assert(not is_poison());
+    return state;
+  }
+  bbox container() const
+  {
+    assert(not is_poison());
+    return bbox();
+  }
+  
+  bool operator==(const bboxset& other) const
+  {
+    assert(not is_poison() and not other.is_poison());
+    return state == other.state;
+  }
+  bool operator<=(const bboxset& other) const
+  {
+    assert(not is_poison() and not other.is_poison());
+    return state <= other.state;
+  }
+  bool operator!=(const bboxset& other) const
+  {
+    return not (*this == other);
+  }
+  bool operator>=(const bboxset& other) const
+  {
+    return other <= *this;
+  }
+  bool operator<(const bboxset& other) const
+  {
+    return not (*this >= other);
+  }
+  bool operator>(const bboxset& other) const
+  {
+    return not (*this <= other);
+  }
+  
+  bboxset& operator&=(const bboxset& other)
+  {
+    assert(not is_poison() and not other.is_poison());
+    state &= other.state;
+    return *this;
+  }
+  bboxset operator&(const bboxset& other) const
+  {
+    bboxset res = *this;
+    return res &= other;
+  }
+  bboxset& operator^=(const bboxset& other)
+  {
+    assert(not is_poison() and not other.is_poison());
+    state ^= other.state;
+    return *this;
+  }
+  bboxset operator^(const bboxset& other) const
+  {
+    bboxset res = *this;
+    return res ^= other;
+  }
+  bboxset& operator|=(const bboxset& other)
+  {
+    assert(not is_poison() and not other.is_poison());
+    state |= other.state;
+    return *this;
+  }
+  bboxset operator|(const bboxset& other) const
+  {
+    bboxset res = *this;
+    return res |= other;
+  }
+  bboxset& operator+=(const bboxset& other)
+  {
+    assert(not is_poison() and not other.is_poison());
+    assert(not (state and other.state));
+    return *this |= other;
+  }
+  bboxset operator+(const bboxset& other) const
+  {
+    bboxset res = *this;
+    return res += other;
+  }
+  bboxset& operator-=(const bboxset& other)
+  {
+    assert(not is_poison() and not other.is_poison());
+    state &= not other.state;
+    return *this;
+  }
+  bboxset operator-(const bboxset& other) const
+  {
+    bboxset res = *this;
+    return res -= other;
+  }
+  
+  bboxset& operator&=(const bbox& b)
+  {
+    return *this &= bboxset(b);
+  }
+  bboxset& operator^=(const bbox& b)
+  {
+    return *this ^= bboxset(b);
+  }
+  bboxset& operator|=(const bbox& b)
+  {
+    return *this |= bboxset(b);
+  }
+  bboxset& operator+=(const bbox& b)
+  {
+    return *this += bboxset(b);
+  }
+  bboxset& operator-=(const bbox& b)
+  {
+    return *this -= bboxset(b);
+  }
+  
+  bboxset shift(const vect& dist, const vect& dist_denom = vect(1)) const
+  {
+    assert(not is_poison());
+    return *this;
+  }
+  bboxset expand(const vect& lo, const vect& hi) const
+  {
+    assert(not is_poison());
+    return *this;
+  }
+  bboxset expand(const ::vect<vect,2>& lohi) const
+  {
+    return expand(lohi[0], lohi[1]);
+  }
+  
+  template<typename C>
+  void serialise(C& out) const
+  {
+    assert(not is_poison());
+    if (state) out.insert(out.end(), bbox());
+  }
+  
+  size_t memory() const
+  {
+    return sizeof *this;
+  }
+  
+  istream& input(istream& is);
+  
+  ostream& output(ostream& os) const;
+};
+
+
+
+/** Operators involving bboxes */
+template<typename T, int D>
+inline bboxset<T,D> operator&(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs & bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator&(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) & bs;
+}
+// Note: bbox & bbox -> bbox
+
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs ^ bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) ^ bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+  return bboxset<T,D>(b1) ^ bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs | bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) | bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+  return bboxset<T,D>(b1) | bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs + bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) + bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+  return bboxset<T,D>(b1) + bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs - bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) - bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+  return bboxset<T,D>(b1) - bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bool operator==(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) == bs;
+}
+template<typename T, int D>
+inline bool operator==(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs == bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator!=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) != bs;
+}
+template<typename T, int D>
+inline bool operator!=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs != bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator<=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) <= bs;
+}
+template<typename T, int D>
+inline bool operator<=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs <= bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator<(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) < bs;
+}
+template<typename T, int D>
+inline bool operator<(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs < bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator>=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) >= bs;
+}
+template<typename T, int D>
+inline bool operator>=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs >= bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator>(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+  return bboxset<T,D>(b) > bs;
+}
+template<typename T, int D>
+inline bool operator>(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+  return bs > bboxset<T,D>(b);
+}
+
+
+
+/** Memory usage */
+template<typename T, int D>
+inline size_t memoryof(bboxset<T,D> const& bs)
+{
+  return bs.memory();
+}
+
+/** Input */
+template<typename T, int D>
+inline istream& operator>>(istream& is, bboxset<T,D>& bs)
+{
+  return bs.input(is);
+}
+
+/** Output */
+template<typename T, int D>
+inline ostream& operator<<(ostream& os, const bboxset<T,D>& bs)
+{
+  return bs.output(os);
+}
+  
+  
+  
+  //////////////////////////////////////////////////////////////////////////////
+  // bboxset<T,D>
+  //////////////////////////////////////////////////////////////////////////////
+  
+  /** Copy constructor */
+  template<typename T, int D>
+  bboxset<T,D>::bboxset(const bboxset& other): 
+    stride(other.stride), offset(other.offset), is_poison_(other.is_poison_)
+  {
+    for (const auto& pos_subset: other.subsets) {
+      const T& pos = pos_subset.first;
+      const auto new_subsetp = make_shared<bboxset1>(*pos_subset.second.get());
+      subsets.insert(subsets.end(), make_pair(pos, new_subsetp));
+    }
+  }
+  
+  /** Copy constructor */
+  // bboxset(bboxset&& other): stride(other.stride), offset(other.offset)
+  // {
+  //   swap(subsets, other.subsets);
+  // }
+  
+  /** Assignment */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator=(const bboxset& other)
+  {
+    if (&other == this) return *this;
+    subsets.clear();
+    for (const auto& pos_subset: other.subsets) {
+      const T& pos = pos_subset.first;
+      const auto new_subsetp = make_shared<bboxset1>(*pos_subset.second.get());
+      subsets.insert(subsets.end(), make_pair(pos, new_subsetp));
+    }
+    stride = other.stride;
+    offset = other.offset;
+    is_poison_ = other.is_poison_;
+    return *this;
+  }
+  
+  /** Assignment */
+  // bboxset& operator=(bboxset&& other)
+  // {
+  //   if (&other == this) return *this;
+  //   swap(subsets, other.subsets);
+  //   stride = other.stride;
+  //   offset = other.offset;
+  //   return *this;
+  // }
+  
+  /** Create set from bbox */
+  template<typename T, int D>
+  bboxset<T,D>::bboxset(const bbox& b):
+    stride(b.stride()), offset(imod(b.lower(), b.stride())), is_poison_(false)
+  {
+    assert(not b.is_poison());
+    if (b.empty()) return;
+    const T lo = last(b.lower());
+    const T hi = last(b.upper()) + last(b.stride());
+    const bbox1 b1(init(b.lower()), init(b.upper()), init(b.stride()));
+    const auto lo_subsetp = make_shared<bboxset1>(b1);
+    const auto hi_subsetp = make_shared<bboxset1>(b1);
+    // subsets.emplace_hint(subsets.end(), lo, lo_subsetp);
+    // subsets.emplace_hint(subsets.end(), hi, hi_subsetp);
+    subsets.insert(subsets.end(), make_pair(lo, lo_subsetp));
+    subsets.insert(subsets.end(), make_pair(hi, hi_subsetp));
+  }
+  
+  /** Create set from container of bboxes, bboxsets, or container
+      thereof */
+  template<typename T, int D>
+  template<typename C>
+  bboxset<T,D>::bboxset(const C& elts)
+  {
+    *this = bboxset();
+    for (const auto& elt: elts) {
+      *this += bboxset(elt);
+    }
+  }
+  
+  /** Create set from container of structs containing a bbox, bboxset,
+      or container thereof */
+  template<typename T, int D>
+  template<typename C, typename S, typename B>
+  bboxset<T,D>::bboxset(const C& elts, const B S::* const mptr)
+  {
+    *this = bboxset();
+    for (const auto& elt: elts) {
+      *this += bboxset(elt.*mptr);
+    }
+  }
+  
+  
+  
+  /** Find a point contained in this set */
+  template<typename T, int D>
+  vect<T,D> bboxset<T,D>::front() const
+  {
+    assert(not is_poison());
+    assert(not empty());
+    const auto& pos_subset = *subsets.begin();
+    const auto& pos = pos_subset.first;
+    const auto& subset = *pos_subset.second.get();
+    const vect1 f1 = subset.front();
+    return vect(f1, pos);
+  }
+  
+  /** Number of characteristic points */
+  template<typename T, int D>
+  T bboxset<T,D>::chi_size() const
+  {
+    assert(not is_poison());
+    T sum = 0;
+    for (const auto& pos_subset: subsets) sum += pos_subset.second->chi_size();
+    return sum;
+  }
+  
+  /** Number of elements */
+  template<typename T, int D>
+  T bboxset<T,D>::size() const
+  {
+    assert(not is_poison());
+    T total_size = 0;           // accumulated total number of points
+    T old_pos = numeric_limits<T>::min(); // location of last subset
+    T old_subset_size = 0;      // number of points in the last subset
+    traverse_subsets
+      ([&](const T& pos, const bboxset1& subset)
+       {
+         const T subset_size = subset.size();
+         total_size += (pos - old_pos) * old_subset_size;
+         old_pos = pos;
+         old_subset_size = subset_size;
+       });
+    assert(old_subset_size == 0);
+    return total_size;
+  }
+  
+  /** Container (min and max) */
+  template<typename T, int D>
+  bbox<T,D> bboxset<T,D>::container() const
+  {
+    assert(not is_poison());
+    if (empty()) return bbox();
+    const T lo = subsets.begin()->first;
+    const T hi = subsets.rbegin()->first;
+    bbox1 container1;
+    traverse_subsets
+      ([&](const T& pos, const bboxset1& subset)
+       {
+         container1 = container1.expanded_containing(subset.container());
+       });
+    return bbox(vect(container1.lower(), lo),
+                vect(container1.upper(), hi - last(stride)),
+                vect(container1.stride(), last(stride)));
+  }
+  
+  
+  
+  /** Test for equality */
+  template<typename T, int D>
+  bool bboxset<T,D>::operator==(const bboxset& other) const
+  {
+    return (*this ^ other).empty();
+  }
+  
+  /** Test for is-subset-of */
+  template<typename T, int D>
+  bool bboxset<T,D>::operator<=(const bboxset& other) const
+  {
+    return (*this | other) == other;
+  }
+  
+  /** Test for is-strict-subset-of */
+  template<typename T, int D>
+  bool bboxset<T,D>::operator<(const bboxset& other) const
+  {
+    return *this != other and *this <= other;
+  }
+  
+  
+  
+  /** Symmetric set difference */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::operator^(const bboxset& other) const
+  {
+    // TODO: If other is much smaller than this, direct insertion may
+    // be faster
+    return binary_operator
+      ([](const bboxset1& set0, const bboxset1& set1) { return set0 ^ set1; },
+       other);
+  }
+  
+  /** Symmetric set difference */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator^=(const bboxset& other)
+  {
+    bboxset res = *this;
+    return *this = *this ^ other;
+  }
+  
+  /** Set intersection */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::operator&(const bboxset& other) const
+  {
+    return binary_operator
+      ([](const bboxset1& set0, const bboxset1& set1) { return set0 & set1; },
+       other);
+  }
+  
+  /** Set intersection */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator&=(const bboxset& other)
+  {
+    return *this = *this & other;
+  }
+  
+  /** Set Union */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::operator|(const bboxset& other) const
+  {
+    return binary_operator
+      ([](const bboxset1& set0, const bboxset1& set1) { return set0 | set1; },
+       other);
+  }
+  
+  /** Set union */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator|=(const bboxset& other)
+  {
+    return *this = *this | other;
+  }
+  
+  /** Symmetric set union */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::operator+(const bboxset& other) const
+  {
+    // return binary_operator
+    //   ([](const bboxset1& set0, const bboxset1& set1) { return set0 + set1; },
+    //    other);
+#ifdef CARPET_DEBUG
+    assert((*this & other).empty());
+#endif
+    return *this ^ other;
+  }
+  
+  /** Symmetric set union */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator+=(const bboxset& other)
+  {
+    return *this = *this + other;
+  }
+  
+  /** Set difference */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::operator-(const bboxset& other) const
+  {
+    return binary_operator
+      ([](const bboxset1& set0, const bboxset1& set1) { return set0 - set1; },
+       other);
+  }
+  
+  /** Set difference */
+  template<typename T, int D>
+  bboxset<T,D>& bboxset<T,D>::operator-=(const bboxset& other)
+  {
+    return *this = *this - other;
+  }
+  
+  
+  
+  /** Shift all points */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::shift(const vect& dist, const vect& dist_denom)
+    const
+  {
+    assert(not is_poison());
+    assert(all(stride % dist_denom == 0));
+    if (all(dist == 0)) return *this;
+    bboxset res;
+    res.stride = stride;
+    res.offset = imod(offset + dist * stride / dist_denom, res.stride);
+    for (const auto& pos_subset: subsets) {
+      const T& pos = pos_subset.first;
+      const bboxset1& subset = *pos_subset.second.get();
+      const T new_pos = pos + last(dist * stride / dist_denom);
+      const auto new_subsetp =
+        make_shared<bboxset1>(subset.shift(init(dist), init(dist_denom)));
+      res.subsets.insert(res.subsets.end(), make_pair(new_pos, new_subsetp));
+    }
+    return res;
+  }
+  
+  /** Expand the set (convolute with a bbox) */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>:: expand(const vect& lo, const vect& hi) const
+  {
+    assert(not is_poison());
+    assert(all(lo >= 0 and hi >= 0));
+    bboxset res = shift(-lo);
+    for (int d=0; d<D; ++d) {
+      T to_expand = (hi + lo)[d];
+      T current_size = 1;
+      while (to_expand > 0) {
+        const T this_expand = min(to_expand, current_size);
+        res |= res.shift(vect::dir(d) * this_expand);
+        current_size += this_expand;
+        to_expand -= this_expand;
+      }
+      assert(to_expand == 0);
+    }
+    return res;
+  }
+  
+  
+  
+  template<typename T, int D>
+  bool bboxset<T,D>::strides_are_compatible(const vect& str1, const vect& str2)
+  {
+    if (all(str1 >= str2)) {
+      return all(imod(str1, str2) == 0);
+    } else if (all(str1 <= str2)) {
+      return all(imod(str2, str1) == 0);
+    }
+    return false;
+  }
+  
+  /** Expande the set (changing the stride): find the smallest
+      b-compatible bboxset around this bboxset */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::expanded_for(const bbox& target) const
+  {
+    // Check preconditions
+    assert(not is_poison() and not target.is_poison());
+    assert(not target.empty());
+    assert(strides_are_compatible(stride, target.stride()));
+    vector<bbox> bs;
+    serialise(bs);
+    bboxset res;
+    for (const auto& b: bs) {
+      res |= b.expanded_for(target);
+    }
+    return res;
+  }
+  
+  /** Shrink the set (changing the stride): find the largest
+      b-compatible bboxset inside this bboxset */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::contracted_for(const bbox& target) const
+  {
+    // Check preconditions
+    assert(not is_poison() and not target.is_poison());
+    assert(not target.empty());
+    assert(strides_are_compatible(stride, target.stride()));
+    const bbox cont = container();
+    const vect safety = 10;
+    const bbox good_world = cont.anti_contracted_for(target).expand(safety);
+    const bbox world1 = good_world.anti_contracted_for(cont).expand(safety);
+    const bbox world2 = world1.anti_contracted_for(target).expand(safety);
+    return (world2 ^ (world1 ^ *this).anti_contracted_for(target)) & good_world;
+  }
+  
+  /** Expande the set (changing the stride): find the smallest open
+      b-compatible bboxset around this bboxset */
+  template<typename T, int D>
+  bboxset<T,D> bboxset<T,D>::anti_contracted_for(const bbox& target) const
+  {
+    // Check preconditions
+    assert(not is_poison() and not target.is_poison());
+    assert(not target.empty());
+    assert(strides_are_compatible(stride, target.stride()));
+    vector<bbox> bs;
+    serialise(bs);
+    bboxset res;
+    for (const auto& b: bs) {
+      res |= b.anti_contracted_for(target);
+    }
+    return res;
+  }
+  
+  
+  
+  /** Serialise the set */
+  template<typename T, int D>
+  template<typename C>
+  void bboxset<T,D>::serialise(C& out) const
+  {
+    assert(not is_poison());
+    typedef map<bbox1,T> subboxes_t;
+    typedef set<bbox1> subboxes1_t;
+    // TODO: Instead of copying from old_subboxes to subboxes,
+    // maintain subboxes via a non-const iterator
+    subboxes_t old_subboxes;
+    traverse_subsets
+      ([&](const T& pos, const bboxset1& subset)
+       {
+         // Convert subset to bboxes
+         subboxes1_t subboxes1;
+         subset.serialise(subboxes1);
+         
+         const bbox1 dummy1;
+         subboxes_t subboxes;
+         // subboxes.reserve(old_subboxes.size() + subboxes1.size());
+         typedef typename subboxes_t::const_iterator subboxes_iter_t;
+         typedef typename subboxes1_t::const_iterator subboxes1_iter_t;
+         subboxes_iter_t iter0 = old_subboxes.begin();
+         subboxes1_iter_t iter1 = subboxes1.begin();
+         subboxes_iter_t const end0 = old_subboxes.end();
+         subboxes1_iter_t const end1 = subboxes1.end();
+         while (iter0!=end0 or iter1!=end1) {
+           bool active0 = iter0!=end0;
+           bool active1 = iter1!=end1;
+           const bbox1& subbox0 = active0 ? iter0->first : dummy1;
+           const bbox1& subbox1 = active1 ? *iter1 : dummy1;
+           // When both subboxes are active, keep only the first (as
+           // determined by less<>)
+           if (active0 and active1) {
+             if (subbox0 != subbox1) {
+               /*const*/ less<bbox1> bbox1_less;
+               if (bbox1_less(subbox0, subbox1)) {
+                 active1 = false;
+               } else {
+                 active0 = false;
+               }
+             }
+           }
+           const T old_pos = active0 ? iter0->second : T();
+           
+           if ((active0 and active1) and (subbox0 == subbox1)) {
+             // The current bbox continues unchanged -- keep it
+             subboxes.insert(subboxes.end(), *iter0);
+           } else {
+             if (active0) {
+               // The current box changed; finalize it
+               const bbox new_box(vect(subbox0.lower(), old_pos),
+                                  vect(subbox0.upper(), pos - last(stride)),
+                                  stride);
+               out.insert(out.end(), new_box);
+             }
+             if (active1) {
+               // There is a new box; add it
+               subboxes.insert(subboxes.end(), make_pair(subbox1, pos));
+             }
+           }
+           
+           if (active0) ++iter0;
+           if (active1) ++iter1;
+         }
+         swap(old_subboxes, subboxes);
+       });
+    assert(old_subboxes.empty());
+  }
+  
+  
+  
+  /** Memory usage */
+  template<typename T, int D>
+  size_t bboxset<T,D>::memory() const
+  {
+    size_t s = sizeof *this;
+    for (const auto& pos_subset: subsets) {
+      s += sizeof pos_subset;
+      const auto* const subsetp = pos_subset.second.get();
+      s += memoryof(*subsetp);
+    }
+    return s;
+  }
+  
+  /** Input */
+  template<typename T, int D>
+  istream& bboxset<T,D>::input(istream& is)
+  {
+    T Tdummy;
+    try {
+      skipws(is);
+      consume(is, "bboxset<");
+      consume(is, typestring(Tdummy));
+      consume(is, ",");
+      int D_;
+      is >> D_;
+      if (D_ != D) {
+        ostringstream msg;
+        msg << "Input error: Wrong bboxset dimension " << D_ << ", expected " << D;
+        CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+        throw input_error();
+      }
+      consume(is, ">(");
+      consume(is, "set<bbox>:");
+      set<bbox> bs;
+      is >> bs;
+      *this = bboxset(bs);
+      consume(is, ",");
+      consume(is, "stride:");
+      is >> stride;
+      consume(is, ",");
+      consume(is, "offset:");
+      is >> offset;
+      consume(is, ")");
+    } catch (input_error& err) {
+      ostringstream msg;
+      msg << "Input error while reading a bboxset<" << typestring(Tdummy) << ",0>";
+      CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+      throw err;
+    }
+    is_poison_ = false;
+    return is;
+  }
+  
+  /** Output */
+  template<typename T, int D>
+  ostream& bboxset<T,D>::output(ostream& os) const
+  {
+    assert(not is_poison());
+    T Tdummy;
+    set<bbox> bs;
+    serialise(bs);
+    return os << "bboxset<" << typestring(Tdummy) << "," << D << ">("
+              << "set<bbox>:" << bs << ","
+              << "stride:" << stride << ","
+              << "offset:" << offset << ")";
+  }
+  
+  
+  
+  //////////////////////////////////////////////////////////////////////////////
+  // bboxset<T,0>
+  //////////////////////////////////////////////////////////////////////////////
+  
+  template<typename T>
+  istream& bboxset<T,0>::input(istream& is)
+  {
+    T Tdummy;
+    try {
+      skipws(is);
+      consume(is, "bboxset<");
+      consume(is, typestring(Tdummy));
+      consume(is, ",0>(");
+      consume(is, "state:");
+      is >> state;
+      consume(is, ",");
+      consume(is, "stride:");
+      is >> stride;
+      consume(is, ",");
+      consume(is, "offset:");
+      is >> offset;
+      consume(is, ")");
+    } catch (input_error& err) {
+      ostringstream msg;
+      msg << "Input error while reading a bboxset<" << typestring(Tdummy) << ",0>";
+      CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+      throw err;
+    }
+    is_poison_ = false;
+    return is;
+  }
+
+
+
+  template<typename T>
+  ostream& bboxset<T,0>::output(ostream& os) const
+  {
+    assert(not is_poison());
+    T Tdummy;
+    return os << "bboxset<" << typestring(Tdummy) << ",0>("
+              << "state:" << state << ","
+              << "stride:" << stride << ","
+              << "offset:" << offset << ")";
+  }
+
+
+
+} // namespace bboxset2
+
+#endif  // #ifndef CARPET_NO_BBOXSET2
+
+
+
+#endif  // #ifndef BBOXSET2_HH
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index d9501aa1f..f2a304e7a 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -6,9 +6,11 @@
 #include <iostream>
 #include <list>
 #include <map>
+#include <memory>
 #include <set>
 #include <stack>
 #include <vector>
+#include <utility>
 
 #include "bbox.hh"
 #include "defs.hh"
@@ -283,6 +285,14 @@ ostream& output (ostream& os, const set<T>& s) {
   return os;
 }
 
+#ifndef CARPET_NO_BBOXSET2
+// Shared pointer output
+template<class T>
+ostream& output (ostream& os, const shared_ptr<T>& s) {
+  return os << "(&" << *s.get() << ")";
+}
+#endif
+
 // Stack output
 template<class T>
 ostream& output (ostream& os, const stack<T>& s) {
@@ -327,13 +337,11 @@ ostream& output (ostream& os, const vector<T>& v) {
 
 template int ipow (int x, int y);
 template CCTK_REAL ipow (CCTK_REAL x, int y);
-template vect<int,dim> ipow (vect<int,dim> x, int y);
+//template vect<int,dim> ipow (vect<int,dim> x, int y);
 template vect<CCTK_REAL,dim> ipow (vect<CCTK_REAL,dim> x, int y);
 
-template bool equals (vector<ibset> const& v, vector<ibset> const& w);
-
-template size_t memoryof (list<ibbox> const & l);
-template size_t memoryof (list<ivect> const & l);
+//template size_t memoryof (list<ibbox> const & l);
+//template size_t memoryof (list<ivect> const & l);
 template size_t memoryof (list<dh*> const & l);
 template size_t memoryof (list<gh*> const & l);
 template size_t memoryof (list<gdata*> const & l);
@@ -342,49 +350,61 @@ template size_t memoryof (list<th*> const & l);
 template size_t memoryof (stack<void*> const & s);
 template size_t memoryof (vector<bool> const & v);
 template size_t memoryof (vector<int> const & v);
-template size_t memoryof (vector<CCTK_REAL> const & v);
-template size_t memoryof (vector<ibbox> const & v);
-template size_t memoryof (vector<ibset> const & v);
+//template size_t memoryof (vector<CCTK_REAL> const & v);
+template size_t memoryof (vector<bbox<int,1> > const & v);
+template size_t memoryof (vector<bbox<int,2> > const & v);
+template size_t memoryof (vector<bbox<int,3> > const & v);
+template size_t memoryof (vector<bboxset1::bboxset<int,dim> > const & v);
+#ifndef CARPET_NO_BBOXSET2
+template size_t memoryof (vector<bboxset2::bboxset<int,dim>> const & v);
+#endif
 template size_t memoryof (vector<ivect> const & v);
 template size_t memoryof (vector<i2vect> const & v);
-template size_t memoryof (vector<fulltree <int,dim,pseudoregion_t> *> const & f);
-template size_t memoryof (vector<pseudoregion_t> const & v);
-template size_t memoryof (vector<region_t> const & v);
+template size_t memoryof (vector<fulltree<int,dim,pseudoregion_t>*> const & f);
+//template size_t memoryof (vector<pseudoregion_t> const & v);
+//template size_t memoryof (vector<region_t> const & v);
 template size_t memoryof (vector<sendrecv_pseudoregion_t> const & v);
 template size_t memoryof (vector<vector<int> > const & v);
 template size_t memoryof (vector<vector<CCTK_REAL> > const & v);
 template size_t memoryof (vector<vector<ibbox> > const & v);
 template size_t memoryof (vector<vector<dh::fast_dboxes> > const & v);
-template size_t memoryof (vector<vector<dh::full_dboxes> > const & v);
+//template size_t memoryof (vector<vector<dh::full_dboxes> > const & v);
 template size_t memoryof (vector<vector<dh::level_dboxes> > const & v);
-template size_t memoryof (vector<vector<dh::light_dboxes> > const & v);
-template size_t memoryof (vector<vector<dh::local_dboxes> > const & v);
-template size_t memoryof (vector<vector<region_t> > const & v);
+//template size_t memoryof (vector<vector<dh::light_dboxes> > const & v);
+//template size_t memoryof (vector<vector<dh::local_dboxes> > const & v);
+//template size_t memoryof (vector<vector<region_t> > const & v);
 template size_t memoryof (vector<vector<vector<CCTK_REAL> > > const & v);
-template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v);
-template size_t memoryof (vector<vector<vector<dh::full_dboxes> > > const & v);
-template size_t memoryof (vector<vector<vector<dh::level_dboxes> > > const & v);
+//template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v);
+//template size_t memoryof (vector<vector<vector<dh::full_dboxes> > > const & v);
+//template size_t memoryof (vector<vector<vector<dh::level_dboxes> > > const & v);
 template size_t memoryof (vector<vector<vector<dh::light_dboxes> > > const & v);
 template size_t memoryof (vector<vector<vector<dh::local_dboxes> > > const & v);
 template size_t memoryof (vector<vector<vector<region_t> > > const & v);
-template size_t memoryof (vector<vector<vector<gdata*> > > const & v);
+//template size_t memoryof (vector<vector<vector<gdata*> > > const & v);
 template size_t memoryof (vector<vector<vector<vector<gdata*> > > > const & v);
 
-template istream& input (istream& os, list<ibbox>& l);
-template istream& input (istream& os, set<ibbox>& s);
+//template istream& input (istream& os, list<ibbox>& l);
+template istream& input (istream& os, set<bbox<int,1> >& s);
+template istream& input (istream& os, set<bbox<int,2> >& s);
+template istream& input (istream& os, set<bbox<int,3> >& s);
 template istream& input (istream& os, vector<int>& v);
-template istream& input (istream& os, vector<CCTK_REAL>& v);
-template istream& input (istream& os, vector<ibbox>& v);
-template istream& input (istream& os, vector<rbbox>& v);
-template istream& input (istream& os, vector<ibset>& v);
+//template istream& input (istream& os, vector<CCTK_REAL>& v);
+template istream& input (istream& os, vector<bbox<int,1> >& v);
+template istream& input (istream& os, vector<bbox<int,2> >& v);
+template istream& input (istream& os, vector<bbox<int,3> >& v);
+//template istream& input (istream& os, vector<rbbox>& v);
+template istream& input (istream& os, vector<bboxset1::bboxset<int,dim> >& v);
+#ifndef CARPET_NO_BBOXSET2
+template istream& input (istream& os, vector<bboxset2::bboxset<int,dim>>& v);
+#endif
 template istream& input (istream& os, vector<ivect>& v);
 template istream& input (istream& os, vector<bbvect>& v);
 template istream& input (istream& os, vector<i2vect>& v);
-template istream& input (istream& os, vector<region_t>& v);
-template istream& input (istream& os, vector<pseudoregion_t>& v);
+//template istream& input (istream& os, vector<region_t>& v);
+//template istream& input (istream& os, vector<pseudoregion_t>& v);
 template istream& input (istream& os, vector<sendrecv_pseudoregion_t>& v);
 template istream& input (istream& os, vector<vector<int> >& v);
-template istream& input (istream& os, vector<vector<CCTK_REAL> >& v);
+//template istream& input (istream& os, vector<vector<CCTK_REAL> >& v);
 template istream& input (istream& os, vector<vector<ibbox> >& v);
 template istream& input (istream& os, vector<vector<rbbox> >& v);
 template istream& input (istream& os, vector<vector<bbvect> >& v);
@@ -393,31 +413,53 @@ template istream& input (istream& os, vector<vector<region_t> >& v);
 template istream& input (istream& os, vector<vector<vector<CCTK_REAL> > >& v);
 template istream& input (istream& os, vector<vector<vector<region_t> > >& v);
 
-template ostream& output (ostream& os, const list<ibbox>& l);
-template ostream& output (ostream& os, const list<region_t>& l);
-template ostream& output (ostream& os, const map<string,Carpet::Timer*>& m);
-template ostream& output (ostream& os, const set<ibbox>& s);
-template ostream& output (ostream& os, const set<ibset>& s);
-template ostream& output (ostream& os, const stack<ibbox>& s);
+//template ostream& output (ostream& os, const list<ibbox>& l);
+//template ostream& output (ostream& os, const list<region_t>& l);
+#ifndef CARPET_NO_BBOXSET2
+//template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,0> >>& m);
+//template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,1> >>& m);
+//template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,2> >>& m);
+//template ostream& output (ostream& os, const map<int,shared_ptr<bboxset2::bboxset<int,3> >>& m);
+#endif
+//template ostream& output (ostream& os, const map<string,Carpet::Timer*>& m);
+//template ostream& output (ostream& os, const pair<bbox<int,0>,int>& p);
+//template ostream& output (ostream& os, const set<bbox<int,0> >& s);
+template ostream& output (ostream& os, const set<bbox<int,1> >& s);
+template ostream& output (ostream& os, const set<bbox<int,2> >& s);
+template ostream& output (ostream& os, const set<bbox<int,3> >& s);
+//template ostream& output (ostream& os, const set<bboxset1::bboxset<int,dim> >& s);
+#ifndef CARPET_NO_BBOXSET2
+//template ostream& output (ostream& os, const shared_ptr<bboxset2::bboxset<int,0> >& s);
+#endif
+//template ostream& output (ostream& os, const stack<ibbox>& s);
 template ostream& output (ostream& os, const vector<bool>& v);
 template ostream& output (ostream& os, const vector<int>& v);
 template ostream& output (ostream& os, const vector<CCTK_REAL>& v);
-template ostream& output (ostream& os, const vector<ibbox>& v);
+template ostream& output (ostream& os, const vector<bbox<int,1> >& v);
+template ostream& output (ostream& os, const vector<bbox<int,2> >& v);
+template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
 template ostream& output (ostream& os, const vector<rbbox>& v);
-template ostream& output (ostream& os, const vector<ibset>& v);
+template ostream& output (ostream& os, const vector<bboxset1::bboxset<int,dim> >& v);
+#ifndef CARPET_NO_BBOXSET2
+template ostream& output (ostream& os, const vector<bboxset2::bboxset<int,dim>>& v);
+#endif
 template ostream& output (ostream& os, const vector<ivect>& v);
 template ostream& output (ostream& os, const vector<rvect>& v);
-template ostream& output (ostream& os, const vector<bbvect>& v);
+//template ostream& output (ostream& os, const vector<bbvect>& v);
 template ostream& output (ostream& os, const vector<i2vect>& v);
-template ostream& output (ostream& os, const vector<dh::fast_dboxes> & v);
-template ostream& output (ostream& os, const vector<dh::full_dboxes> & v);
-template ostream& output (ostream& os, const vector<dh::level_dboxes> & v);
-template ostream& output (ostream& os, const vector<dh::light_dboxes> & v);
-template ostream& output (ostream& os, const vector<dh::local_dboxes> & v);
+//template ostream& output (ostream& os, const vector<dh::fast_dboxes> & v);
+//template ostream& output (ostream& os, const vector<dh::full_dboxes> & v);
+//template ostream& output (ostream& os, const vector<dh::level_dboxes> & v);
+//template ostream& output (ostream& os, const vector<dh::light_dboxes> & v);
+//template ostream& output (ostream& os, const vector<dh::local_dboxes> & v);
 template ostream& output (ostream& os, const vector<region_t>& v);
-template ostream& output (ostream& os, const vector<pseudoregion_t>& v);
+//template ostream& output (ostream& os, const vector<pseudoregion_t>& v);
 template ostream& output (ostream& os, const vector<sendrecv_pseudoregion_t>& v);
-template ostream& output (ostream& os, const vector<list<ibbox> >& v);
+//template ostream& output (ostream& os, const vector<list<ibbox> >& v);
+//template ostream& output (ostream& os, const vector<pair<bbox<int,0>,int> >& v);
+//template ostream& output (ostream& os, const vector<pair<bbox<int,1>,int> >& v);
+//template ostream& output (ostream& os, const vector<pair<bbox<int,2>,int> >& v);
+//template ostream& output (ostream& os, const vector<pair<bbox<int,3>,int> >& v);
 template ostream& output (ostream& os, const vector<vector<int> >& v);
 template ostream& output (ostream& os, const vector<vector<CCTK_REAL> >& v);
 template ostream& output (ostream& os, const vector<vector<ibbox> >& v);
@@ -425,16 +467,16 @@ template ostream& output (ostream& os, const vector<vector<rbbox> >& v);
 template ostream& output (ostream& os, const vector<vector<bbvect> >& v);
 template ostream& output (ostream& os, const vector<vector<i2vect> >& v);
 template ostream& output (ostream& os, const vector<vector<dh::fast_dboxes> > & b);
-template ostream& output (ostream& os, const vector<vector<dh::full_dboxes> > & b);
+//template ostream& output (ostream& os, const vector<vector<dh::full_dboxes> > & b);
 template ostream& output (ostream& os, const vector<vector<dh::level_dboxes> > & b);
-template ostream& output (ostream& os, const vector<vector<dh::light_dboxes> > & b);
-template ostream& output (ostream& os, const vector<vector<dh::local_dboxes> > & b);
+//template ostream& output (ostream& os, const vector<vector<dh::light_dboxes> > & b);
+//template ostream& output (ostream& os, const vector<vector<dh::local_dboxes> > & b);
 template ostream& output (ostream& os, const vector<vector<region_t> >& v);
 template ostream& output (ostream& os, const vector<vector<vector<CCTK_REAL> > >& v);
-template ostream& output (ostream& os, const vector<vector<vector<ibbox> > >& v);
-template ostream& output (ostream& os, const vector<vector<vector<dh::fast_dboxes> > > & b);
-template ostream& output (ostream& os, const vector<vector<vector<dh::full_dboxes> > > & b);
-template ostream& output (ostream& os, const vector<vector<vector<dh::level_dboxes> > > & b);
+//template ostream& output (ostream& os, const vector<vector<vector<ibbox> > >& v);
+//template ostream& output (ostream& os, const vector<vector<vector<dh::fast_dboxes> > > & b);
+//template ostream& output (ostream& os, const vector<vector<vector<dh::full_dboxes> > > & b);
+//template ostream& output (ostream& os, const vector<vector<vector<dh::level_dboxes> > > & b);
 template ostream& output (ostream& os, const vector<vector<vector<dh::light_dboxes> > > & b);
 template ostream& output (ostream& os, const vector<vector<vector<dh::local_dboxes> > > & b);
 template ostream& output (ostream& os, const vector<vector<vector<region_t> > >& v);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 7968300a2..e52c94d12 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -10,6 +10,7 @@
 #include <iostream>
 #include <list>
 #include <map>
+#include <memory>
 #include <set>
 #include <stack>
 #include <vector>
@@ -62,13 +63,16 @@ char const * const eol = "\n";
 #else
 #  define AT(index) operator[](index)
 #endif
-  
+
 
 
 // Some shortcuts for type names
 template<typename T, int D> class vect;
 template<typename T, int D> class bbox;
-template<typename T, int D> class bboxset;
+namespace bboxset1 { template<typename T, int D> class bboxset; }
+#ifndef CARPET_NO_BBOXSET2
+namespace bboxset2 { template<typename T, int D> class bboxset; }
+#endif
 template<typename T, int D, typename P> class fulltree;
 
 typedef vect<bool,dim>      bvect;
@@ -78,15 +82,18 @@ typedef vect<CCTK_REAL,dim> rvect;
 typedef bbox<int,dim>       ibbox;
 typedef bbox<CCTK_INT,dim>  jbbox;
 typedef bbox<CCTK_REAL,dim> rbbox;
-typedef bboxset<int,dim>    ibset;
-  
+namespace bboxset1 { typedef bboxset<int,dim> ibset; }
+#ifndef CARPET_NO_BBOXSET2
+namespace bboxset2 { typedef bboxset<int,dim> ibset; }
+#endif
+
 // (Try to replace these by b2vect and i2vect)
-typedef vect<vect<bool,2>,dim> bbvect;
-typedef vect<vect<int,2>,dim>  iivect;
+typedef vect<vect<bool,2>,dim>     bbvect;
+typedef vect<vect<int,2>,dim>      iivect;
 typedef vect<vect<CCTK_INT,2>,dim> jjvect;
 
-typedef vect<vect<bool,dim>,2> b2vect;
-typedef vect<vect<int,dim>,2>  i2vect;
+typedef vect<vect<bool,dim>,2>      b2vect;
+typedef vect<vect<int,dim>,2>       i2vect;
 typedef vect<vect<CCTK_INT,dim>,2>  j2vect;
 typedef vect<vect<CCTK_REAL,dim>,2> r2vect;
 
@@ -331,6 +338,9 @@ template<class T> ostream& output (ostream& os, const list<T>& l);
 template<class S, class T> ostream& output (ostream& os, const map<S,T>& m);
 template<class S, class T> ostream& output (ostream& os, const pair<S,T>& p);
 template<class T> ostream& output (ostream& os, const set<T>& s);
+#ifndef CARPET_NO_BBOXSET2
+template<class T> ostream& output (ostream& os, const shared_ptr<T>& s);
+#endif
 template<class T> ostream& output (ostream& os, const stack<T>& s);
 template<class T> ostream& output (ostream& os, const vector<T>& v);
 
@@ -344,11 +354,23 @@ inline ostream& operator<< (ostream& os, const map<S,T>& m) {
   return output(os,m);
 }
 
+template<class S, class T>
+inline ostream& operator<< (ostream& os, const pair<S,T>& s) {
+  return output(os,s);
+}
+
 template<class T>
 inline ostream& operator<< (ostream& os, const set<T>& s) {
   return output(os,s);
 }
 
+#ifndef CARPET_NO_BBOXSET2
+template<class T>
+inline ostream& operator<< (ostream& os, const shared_ptr<T>& s) {
+  return output(os,s);
+}
+#endif
+
 template<class T>
 inline ostream& operator<< (ostream& os, const stack<T>& s) {
   return output(os,s);
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 01a4210bb..7d5afb7d8 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -5,6 +5,8 @@ SRCS =	backtrace.cc				\
 	balance.cc				\
 	bbox.cc					\
 	bboxset.cc				\
+	bboxset1.cc				\
+	bboxset2.cc				\
 	bintree.cc				\
 	cacheinfo.cc				\
 	commstate.cc				\
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
index f28fe2d81..cfc65e0ba 100644
--- a/Carpet/CarpetLib/src/region.hh
+++ b/Carpet/CarpetLib/src/region.hh
@@ -7,6 +7,7 @@
 #include "defs.hh"
 #include "dist.hh"
 #include "bbox.hh"
+#include "bboxset.hh"
 #include "fulltree.hh"
 #include "vect.hh"