Introduce a tree data structure to speed up domain decomposition

Introduce a tree data structure "fulltree", which decomposes a single,
rectangular region into a tree of non-overlapping, rectangular sub-regions.

Move the processor decomposition from the regridding thorns into Carpet.
Create such trees during processor decomposition.
Store these trees with the grid hierarchy.

12 files changed, 925 insertions, 33 deletions

diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index 98a160bf9..44b7729d5 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -7,6 +7,7 @@ includes header: dist.hh in dist.hh
 
 includes header: bbox.hh in bbox.hh
 includes header: bboxset.hh in bboxset.hh
+includes header: fulltree.hh in fulltree.hh
 includes header: region.hh in region.hh
 includes header: vect.hh in vect.hh
 
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 13ef01301..357061ee2 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -224,6 +224,7 @@ ostream& output (ostream& os, const vector<T>& v) {
 #include "bbox.hh"
 #include "bboxset.hh"
 #include "dh.hh"
+#include "fulltree.hh"
 #include "gdata.hh"
 #include "ggf.hh"
 #include "region.hh"
@@ -245,6 +246,7 @@ template size_t memoryof (vector<int> const & v);
 template size_t memoryof (vector<CCTK_REAL> const & v);
 template size_t memoryof (vector<bbox<int,3> > const & v);
 template size_t memoryof (vector<vect<int,3> > const & v);
+template size_t memoryof (vector<fulltree <int,3,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);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 59d22f4c8..1152c154a 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -57,14 +57,19 @@ const int dim = 3;
 
 
 // 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;
-template<typename T, int D> class vect;
+template<typename T, int D, typename P> class fulltree;
+
+struct pseudoregion_t;
+struct region_t;
 
 typedef vect<bool,dim>   bvect;
 typedef vect<int,dim>    ivect;
 typedef bbox<int,dim>    ibbox;
 typedef bboxset<int,dim> ibset;
+typedef fulltree<int,dim,pseudoregion_t> ipfulltree;
 
 // (Try to replace these by b2vect and i2vect)
 typedef vect<vect<bool,2>,dim> bbvect;
diff --git a/Carpet/CarpetLib/src/fulltree.cc b/Carpet/CarpetLib/src/fulltree.cc
new file mode 100644
index 000000000..b114dee19
--- /dev/null
+++ b/Carpet/CarpetLib/src/fulltree.cc
@@ -0,0 +1,492 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+
+#include "defs.hh"
+#include "region.hh"
+
+#include "fulltree.hh"
+
+
+
+#if 0
+// Create an empty tree
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree ()
+  : type (type_empty)
+{
+  assert (invariant());
+}
+#endif
+
+
+
+// Create a tree branch from a list of bounds and subtrees
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (int dir_, vector <T> const & bounds_,
+                           vector <fulltree *> const & subtrees_)
+  : type (type_branch), dir (dir_), bounds (bounds_), subtrees (subtrees_)
+{
+  assert (dir>=0 and dir<D);
+  assert (bounds.size() == subtrees.size() + 1);
+  assert (invariant());
+}
+
+
+
+// Create a tree leaf from a payload
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (P const & p_)
+  : type (type_leaf), p (p_)
+{
+  assert (invariant());
+}
+  
+
+
+// Create a tree as copy from another tree
+template <typename T, int D, typename P>
+fulltree<T,D,P>::fulltree (fulltree const & t)
+  : type (t.type)
+{
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    subtrees.resize (t.subtrees.size());
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      subtrees.AT(i) = new fulltree (*t.subtrees.AT(i));
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+}
+  
+
+
+// Copy a tree
+template <typename T, int D, typename P>
+fulltree<T,D,P> &
+fulltree<T,D,P>::operator= (fulltree const & t)
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      delete subtrees.AT(i);
+    }
+  }
+  type = t.type;
+  switch (type) {
+  case type_empty:
+    // do nothing
+    break;
+  case type_branch:
+    dir = t.dir;
+    bounds = t.bounds;
+    subtrees.resize (t.subtrees.size());
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      subtrees.AT(i) = new fulltree (*t.subtrees.AT(i));
+    }
+    break;
+  case type_leaf:
+    p = t.p;
+    break;
+  default:
+    assert (0);
+  }
+  assert (invariant());
+  return *this;
+}
+
+
+
+// Delete a tree (and its subtrees)
+template <typename T, int D, typename P>
+fulltree<T,D,P>::~fulltree ()
+{
+  assert (invariant());
+  if (is_branch()) {
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      delete subtrees.AT(i);
+    }
+  }
+}
+
+
+
+// Compare trees
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::operator== (fulltree const & t) const
+{
+  assert (invariant());
+  assert (t.invariant());
+  if (type != t.type) return false;
+  switch (type) {
+  case type_empty:
+    break;
+  case type_branch:
+    if (dir != t.dir) return false;
+    if (bounds != t.bounds) return false;
+    if (subtrees != t.subtrees) return false;
+    break;
+  case type_leaf:
+    return p == t.p;
+  default:
+    assert (0);
+  }
+  return true;
+}
+
+
+
+// Invariant
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::invariant () const
+{
+  return empty() or is_branch() or is_leaf();
+}
+
+
+
+// Find the leaf payload corresponding to a position
+template <typename T, int D, typename P>
+P const *
+fulltree<T,D,P>::search (tvect const & ipos) const
+{
+  assert (not empty());
+  // cout << "fulltree::search ipos=" << ipos << endl;
+  if (is_leaf()) return & p;
+  int const i = asearch (ipos[dir], bounds);
+  // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl;
+  if (i<0 or i>=int(subtrees.size())) return NULL; // not found
+  return subtrees.AT(i)->search(ipos);
+}
+
+template <typename T, int D, typename P>
+P *
+fulltree<T,D,P>::search (tvect const & ipos)
+{
+  assert (not empty());
+  // cout << "fulltree::search ipos=" << ipos << endl;
+  if (is_leaf()) return & p;
+  int const i = asearch (ipos[dir], bounds);
+  // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl;
+  if (i<0 or i>=int(subtrees.size())) return NULL; // not found
+  return subtrees.AT(i)->search(ipos);
+}
+
+
+
+// Const iterator
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    assert (f.subtrees.size() > 0);
+    it = new const_iterator (* f.subtrees.at(i));
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = f.subtrees.size();
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::const_iterator::~const_iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P> const & 
+fulltree<T,D,P>::const_iterator::operator* () const
+{
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::const_iterator::operator== (const_iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename fulltree<T,D,P>::const_iterator &
+fulltree<T,D,P>::const_iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    if ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (not done()) {
+        // to do: use a new function "reset iterator" instead
+        it = new const_iterator (* f.subtrees.at(i));
+      }
+    }
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::const_iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == f.subtrees.size();
+  }
+}
+
+
+
+// Non-const iterator
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::iterator (fulltree & f_)
+  : f(f_), i(0), it(0)
+{
+  if (f.is_branch()) {
+    assert (f.subtrees.size() > 0);
+    it = new iterator (* f.subtrees.at(i));
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::iterator (fulltree & f_, int)
+  : f(f_), it(0)
+{
+  if (f.empty()) {
+    i = 0;
+  } else if (f.is_leaf()) {
+    i = 1;
+  } else {
+    i = f.subtrees.size();
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P>::iterator::~iterator ()
+{
+  if (it) {
+    delete it;
+  }
+}
+
+template <typename T, int D, typename P>
+fulltree<T,D,P> & 
+fulltree<T,D,P>::iterator::operator* ()
+{
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    return f;
+  } else {
+    return **it;
+  }
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::iterator::operator== (iterator const & it2) const
+{
+  // assert (f == it2.f);
+  assert (&f == &it2.f);
+  if (i != it2.i) return false;
+  if (it == 0 and it2.it == 0) return true;
+  if (it == 0 or it2.it == 0) return false;
+  return *it == *it2.it;
+}
+
+template <typename T, int D, typename P>
+typename fulltree<T,D,P>::iterator &
+fulltree<T,D,P>::iterator::operator++ ()
+{
+  assert (not done());
+  assert (not f.empty());
+  if (f.is_leaf()) {
+    ++ i;
+  } else {
+    ++ *it;
+    if ((*it).done()) {
+      delete it;
+      it = 0;
+      ++ i;
+      if (not done()) {
+        // to do: use a new function "reset iterator" instead
+        it = new iterator (* f.subtrees.at(i));
+      }
+    }
+  }
+  return *this;
+}
+
+template <typename T, int D, typename P>
+bool
+fulltree<T,D,P>::iterator::done () const
+{
+  if (f.empty()) {
+    return true;
+  } else if (f.is_leaf()) {
+    return i > 0;
+  } else {
+    return i == f.subtrees.size();
+  }
+}
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+size_t
+fulltree<T,D,P>::memory () const
+{
+  size_t size = sizeof *this;
+  if (is_branch()) {
+    size += memoryof (bounds) + memoryof (subtrees);
+    for (typename vector <fulltree *>::const_iterator
+           i = subtrees.begin(); i != subtrees.end(); ++ i)
+    {
+      size += memoryof(**i);
+    }
+  }
+  return size;
+}
+
+
+
+// Output helper
+template <typename T, int D, typename P>
+void
+fulltree<T,D,P>::output (ostream & os) const
+{
+  os << "fulltree{";
+  if (empty()) {
+    os << "empty";
+  } else if (is_branch()) {
+    os << "branch:"
+       << "dir=" << dir << ","
+       << "subtrees=[";
+    for (size_t i=0; i<subtrees.size(); ++i) {
+      os << bounds.at(i) << ":[" << i << "]=" << subtrees.at(i) << ":";
+    }
+    os << bounds.at(subtrees.size()) << "]";
+  } else {
+    os << "leaf:"
+       << "payload=" << p;
+  }
+  os << "}";
+}
+
+
+
+// Generic arithmetic search
+template <typename T>
+static
+int asearch (T const t, vector <T> const & ts)
+{
+  // cout << "fulltree::asearch t=" << t << " ts=" << ts << endl;
+  int imin = 0;
+  int imax = int(ts.size()) - 1;
+  if (imax <= imin) {
+    // cout << "fulltree::asearch: no values" << endl;
+    return imin;
+  }
+  T tmin = ts.AT(imin);
+  // cout << "fulltree::asearch: imin=" << imin << " tmin=" << tmin << endl;
+  if (t < tmin) {
+    // cout << "fulltree::asearch: below minimum" << endl;
+    return -1;
+  }
+  T tmax = ts.AT(imax);
+  // cout << "fulltree::asearch: imax=" << imax << " tmax=" << tmax << endl;
+  if (t >= tmax) {
+    // cout << "fulltree::asearch: above maximum" << endl;
+    return imax;
+  }
+  int isize = imax - imin;
+  for (;;) {
+    // check loop invariants
+    assert (imin < imax);
+    assert (t >= tmin);
+    assert (t < tmax);
+    // cout << "fulltree::asearch t=" << t << " imin=" << imin << " imax=" << imax << endl;
+    if (imax == imin+1) {
+      // cout << "fulltree::asearch: found value" << endl;
+      return imin;
+    }
+    assert (tmax > tmin);       // require that ts is strictly ordered
+    CCTK_REAL const rguess =
+      (imax - imin) * CCTK_REAL(t - tmin) / (tmax - tmin);
+    int const iguess = imin + max (1, int(floor(rguess)));
+    // handle round-off errors
+    if (iguess == imax) {
+      // cout << "fulltree::asearch: accidental hit after roundoff error" << endl;
+      return imax - 1;
+    }
+    assert (iguess > imin and iguess < imax);
+    T const tguess = ts.AT(iguess);
+    // cout << "fulltree::asearch: intersecting at iguess=" << iguess << " tguess=" << tguess << endl;
+    if (t < tguess) {
+      imax = iguess;
+      tmax = tguess;
+      // cout << "fulltree::asearch: new imax=" << imax << " tmax=" << tmax << endl;
+    } else {
+      imin = iguess;
+      tmin = tguess;
+      // cout << "fulltree::asearch: new imin=" << imin << " tmin=" << tmin << endl;
+    }
+    // check loop variant
+    int const newisize = imax - imin;
+    assert (newisize < isize);
+    isize = newisize;
+  }
+}
+
+
+
+template class fulltree <int, dim, pseudoregion_t>;
+
+template size_t memoryof (fulltree <int, dim, pseudoregion_t> const & f);
+template size_t memoryof (fulltree <int, dim, pseudoregion_t> * const & f);
+
+template ostream & operator<< (ostream & os, fulltree <int, dim, pseudoregion_t> const & f);
diff --git a/Carpet/CarpetLib/src/fulltree.hh b/Carpet/CarpetLib/src/fulltree.hh
new file mode 100644
index 000000000..6ffca209e
--- /dev/null
+++ b/Carpet/CarpetLib/src/fulltree.hh
@@ -0,0 +1,181 @@
+#ifndef FULLTREE_HH
+#define FULLTREE_HH
+
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+#include <vect.hh>
+
+using namespace std;
+
+
+
+// This is a "full tree" data structure, i.e., a tree data structure
+// which decomposes a cuboid domain into a set of non-overlapping
+// cuboid subdomains.  It is an n-ary tree, i.e., each tree node can
+// have arbitrarily many subtrees.  Each node splits a domain in
+// exactly one direction.  Subdomains cannot be empty.
+
+// All intervals are closed at the lower and open at the upper
+// boundary.  This makes it possible to combine adjacent such
+// intervals, obtaining again an interval with this property.  (In
+// particular, if bboxes are stored, the upper bound of bboxes must be
+// increased by the stride to arrive at such intervals.)
+
+
+// Generic arithmetic search
+template <typename T>
+static
+int asearch (T t, vector <T> const & ts);
+
+
+
+// T: index space (usually integer, or maybe real)
+// D: number of dimensions (rank)
+// P: payload (e.g. region_t)
+template <typename T, int D, typename P>
+class fulltree {
+  
+public:
+  
+  // Short name for a small vector
+  typedef vect <T, D> tvect;
+  
+private:
+  
+  enum tree_type_t { type_empty, type_branch, type_leaf } type;
+  
+  // Direction in which the node is split (0 <= dir < D)
+  int dir;
+  
+  // If this is a branch:
+  // n+1 bounds, splitting the domain
+  vector <T> bounds;            // [n+1]
+  // n pointers to subtrees
+  vector <fulltree *> subtrees; // [n]
+  
+  // If this is a leaf:
+  // just the payload
+  P p;
+  
+public:
+  
+#if 0
+  // Create an empty tree
+  fulltree ();
+#endif
+  
+  // Create a tree branch from a list of bounds and subtrees
+  fulltree (int dir_, vector <T> const & bounds_,
+            vector <fulltree *> const & subtrees_);
+  
+  // Create a tree leaf from a payload
+  fulltree (P const & p_);
+  
+  // Create a tree as copy from another tree
+  fulltree (fulltree const & t);
+  
+  // Copy a tree
+  fulltree & operator= (fulltree const & t);
+  
+  // Delete a tree (and its subtrees)
+  ~fulltree ();
+  
+  // Inquire tree properties
+  bool empty() const { return type == type_empty; }
+  bool is_branch() const { return type == type_branch; }
+  bool is_leaf() const { return type == type_leaf; }
+  
+  // Compare trees
+  bool operator== (fulltree const & t) const;
+  bool operator!= (fulltree const & t) const
+  { return not (*this == t); }
+  
+  // Invariant
+  bool invariant () const;
+  
+  // Access the payload
+  P const & payload () const { assert (is_leaf()); return p; }
+  P & payload () { assert (is_leaf()); return p; }
+  
+  // Find the leaf payload corresponding to a position
+  P const * search (tvect const & ipos) const;
+  P * search (tvect const & ipos);
+  
+  class const_iterator {
+    fulltree const & f;
+    size_t i;
+    const_iterator * it;
+  public:
+    const_iterator (fulltree const & f_);
+    const_iterator (fulltree const & f_, int);
+    ~const_iterator ();
+    fulltree const & operator* () const;
+    bool operator== (const_iterator const & it2) const;
+    bool operator!= (const_iterator const & it2) const
+    { return not (*this == it2); }
+    const_iterator & operator++ ();
+    bool done () const;
+  };
+  
+  class iterator {
+    fulltree & f;
+    size_t i;
+    iterator * it;
+  public:
+    iterator (fulltree & f_);
+    iterator (fulltree & f_, int);
+    ~iterator ();
+    fulltree & operator* ();
+    bool operator== (iterator const & it2) const;
+    bool operator!= (iterator const & it2) const
+    { return not (*this == it2); }
+    iterator & operator++ ();
+    bool done () const;
+  };
+  
+  const_iterator begin() const
+  { return const_iterator (*this); }
+  
+  const_iterator end() const
+  { return const_iterator (*this, 0); }
+  
+  iterator begin()
+  { return iterator (*this); }
+  
+  iterator end()
+  { return iterator (*this, 0); }
+  
+  // Memory usage
+  size_t memory () const;
+  
+  // Output helper
+  void output (ostream & os) const;
+};
+
+
+
+// Memory usage
+template <typename T, int D, typename P>
+inline size_t memoryof (fulltree<T,D,P> const & f) { return f.memory(); }
+
+template <typename T, int D, typename P>
+inline size_t memoryof (fulltree<T,D,P> * const & fp) { return sizeof fp; }
+
+
+
+// Output
+template <typename T, int D, typename P>
+ostream &
+operator<< (ostream & os, fulltree<T,D,P> const & f)
+{
+  f.output (os);
+  return os;
+}
+
+
+
+#endif  // #ifndef FULLTREE_HH
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index b6eefa499..ae42668d8 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -512,8 +512,8 @@ transfer_from_all (comm_state & state,
     pseudoregion_t const & precv = (* ipsendrecv).recv;
     ibbox const & send = psend.extent;
     ibbox const & recv = precv.extent;
-    int const c2 = psend.processor;
-    int const c1 = precv.processor;
+    int const c2 = psend.component;
+    int const c1 = precv.component;
     
     // Source and destination data
     gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1);
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index e3c351b9d..b2c86b8db 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -62,7 +62,7 @@ protected:
 public:
   const int vectorlength;       // vector length
   const int vectorindex;        // index of *this
-  const ggf* vectorleader;      // first vector element
+  ggf* const vectorleader;      // first vector element
   
 private:
   mdata oldstorage;             // temporary storage
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index fa9ff2125..9b4a5f443 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -15,6 +15,8 @@
 
 using namespace std;
 
+using namespace CarpetLib;
+
 
 
   // Constructors
@@ -76,10 +78,12 @@ gh::
 // Modifiers
 void
 gh::
-regrid (mregs const & regs)
+regrid (rregs const & superregs, mregs const & regs)
 {
   DECLARE_CCTK_PARAMETERS;
   
+  superregions = superregs;
+  
   // Save the grid hierarchy
   oldregions.clear ();
   swap (oldregions, regions);
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 0a283b4cf..ffac38804 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -47,6 +47,10 @@ public:				// should be readonly
   vector<vector<ibbox> > baseextents; // [ml][rl]
   const i2vect boundary_width;
   
+  // Extents of the regions before distributing them over the
+  // processors
+  rregs superregions;
+  
   mregs regions;                // extents and properties of all grids
   mregs oldregions;             // a copy, used during regridding
   
@@ -65,7 +69,7 @@ public:
   ~gh ();
   
   // Modifiers
-  void regrid (mregs const & regs);
+  void regrid (rregs const & superregs, mregs const & regs);
   bool recompose (int rl, bool do_prolongate);
   
 private:
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index aa9b399c0..88f6261ce 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -8,6 +8,7 @@ SRCS =	bbox.cc					\
 	defs.cc					\
 	dh.cc					\
 	dist.cc					\
+	fulltree.cc				\
 	gdata.cc				\
 	gf.cc					\
 	ggf.cc					\
diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc
index caf535598..0230d373d 100644
--- a/Carpet/CarpetLib/src/region.cc
+++ b/Carpet/CarpetLib/src/region.cc
@@ -1,5 +1,8 @@
+#include <cassert>
+#include <cstdlib>
 #include <iostream>
 
+#include "bboxset.hh"
 #include "defs.hh"
 #include "region.hh"
 
@@ -7,6 +10,69 @@ using namespace std;
 
 
 
+region_t::region_t ()
+  : processor (-1), processors (NULL)
+{
+  assert (invariant());
+}
+
+region_t::region_t (region_t const & a)
+{
+  assert (a.invariant());
+  extent = a.extent;
+  outer_boundaries = a.outer_boundaries;
+  map = a.map;
+  processor = a.processor;
+  if (a.processors == NULL) {
+    processors = NULL;
+  } else {
+    processors = new ipfulltree (*a.processors);
+  }
+  assert (invariant());
+}
+
+region_t &
+region_t::operator= (region_t const & a)
+{
+  assert (invariant());
+  if (processors != NULL) {
+    delete processors;
+  }
+  assert (a.invariant());
+  extent = a.extent;
+  outer_boundaries = a.outer_boundaries;
+  map = a.map;
+  processor = a.processor;
+  if (a.processors == NULL) {
+    processors = NULL;
+  } else {
+    processors = new ipfulltree (*a.processors);
+  }
+  assert (invariant());
+  return *this;
+}
+
+
+region_t::~region_t ()
+{
+  assert (invariant());
+  if (processors != NULL) {
+    delete processors;
+  }
+}
+
+
+
+bool
+region_t::invariant () const
+{
+  if (processor >= 0 and processors != NULL) return false;
+  return true;
+}
+
+
+
+// Compare two regions for equality.
 bool
 operator== (region_t const & a, region_t const & b)
 {
@@ -14,7 +80,113 @@ operator== (region_t const & a, region_t const & b)
     a.extent == b.extent and
     all (all (a.outer_boundaries == b.outer_boundaries)) and
     a.map == b.map and
-    a.processor == b.processor;
+    a.processor == b.processor and
+    ((a.processors == NULL and b.processors == NULL) or
+     (a.processors != NULL and b.processors != NULL and
+      *a.processors == *b.processors));
+}
+
+
+
+// Combine a collection of regions.  Regions can be combined if they
+// abutt on boundaries which are not outer boundaries, ignoring the
+// processor distribution.  This should lead to a canonical
+// representations of collections of regions.
+//
+// We use vectors to represent the collection, but we could also use
+// other containers.  oldregs is read, newregs is added-to.  newregs
+// is not cleared.
+void
+combine_regions (vector<region_t> const & oldregs,
+                 vector<region_t> & newregs)
+{
+  // Find the union of all regions' bounding boxes, and the union of
+  // all regions' outer boundaries.  Represent the boundaries as the
+  // outermost layer of grid points of the corresponding bounding
+  // boxes.
+  int const m = oldregs.empty() ? -1 : oldregs.AT(0).map;
+  ibset comps;
+  ibset cobnds[2][dim];
+  for (vector<region_t>::const_iterator
+         ri = oldregs.begin(); ri != oldregs.end(); ++ ri)
+  {
+    region_t const & reg = * ri;
+    assert (reg.map == m);
+    comps += reg.extent;
+    for (int f = 0; f < 2; ++ f) {
+      for (int d = 0; d < dim; ++ d) {
+        if (reg.outer_boundaries[f][d]) {
+          ibbox bnd = reg.extent;
+          ivect lo = bnd.lower();
+          ivect up = bnd.upper();
+          if (f==0) {
+            up[d] = lo[d];
+          } else {
+            lo[d] = up[d];
+          }
+          bnd = ibbox (lo, up, bnd.stride());
+          cobnds[f][d] += bnd;
+        }
+      }
+    }
+  }
+  comps.normalize();
+  for (int f = 0; f < 2; ++ f) {
+    for (int d = 0; d < dim; ++ d) {
+      cobnds[f][d].normalize();
+    }
+  }
+  
+  // Reserve (generous) memory for the result
+  size_t const needsize = newregs.size() + comps.setsize();
+  if (newregs.capacity() < needsize) {
+    newregs.reserve (1000 + 2 * needsize);
+  }
+  
+  // Insert the regions
+  for (ibset::const_iterator ci = comps.begin(); ci != comps.end(); ++ ci) {
+    ibbox const & c = * ci;
+    b2vect obnds;
+    for (int f = 0; f < 2; ++ f) {
+      for (int d = 0; d < dim; ++ d) {
+        obnds[f][d] = cobnds[f][d].intersects (c);
+        if (obnds[f][d]) {
+          ivect lo = c.lower();
+          ivect up = c.upper();
+          if (f) lo[d]=up[d]; else up[d]=lo[d];
+          ibbox const cbnds (lo, up, c.stride());
+          if (not ((cobnds[f][d] & ibset(c)) == ibset(cbnds))) {
+            cout << "cobnds[f][d] = " << cobnds[f][d] << endl
+                 << "ibset(c) = " << ibset(c) << endl
+                 << "(cobnds[f][d] & ibset(c)) = " << (cobnds[f][d] & ibset(c)) << endl
+                 << "ibset(cbnds) = " << ibset(cbnds) << endl;
+          }
+          assert ((cobnds[f][d] & ibset(c)) == ibset(cbnds));
+        }
+      }
+    }
+    
+    region_t reg;
+    reg.extent           = c;
+    reg.outer_boundaries = obnds;
+    reg.map              = m;
+    reg.processor        = -1;
+    reg.processors       = NULL;
+    newregs.push_back (reg);
+  }
+}
+
+
+
+size_t memoryof (region_t const & reg)
+{
+  return
+    memoryof (reg.extent) +
+    memoryof (reg.outer_boundaries) +
+    memoryof (reg.map) +
+    memoryof (reg.processor) +
+    memoryof (reg.processors) +
+    (reg.processors != NULL ? memoryof (*reg.processors) : 0);
 }
 
 
@@ -59,6 +231,8 @@ operator>> (istream & is, region_t & reg)
   skipws (is);
   consume (is, ')');
   
+  reg.processors = NULL;
+  
   return is;
 }
 
@@ -77,9 +251,20 @@ operator<< (ostream & os, region_t const & reg)
 
 
 
+// Compare two pseudoregions for equality.
+bool
+operator== (pseudoregion_t const & a, pseudoregion_t const & b)
+{
+  return
+    a.extent == b.extent and
+    a.component == b.component;
+}
+
+
+
 ostream & operator<< (ostream & os, pseudoregion_t const & p)
 {
-  return os << p.extent << "/p:" << p.processor;
+  return os << p.extent << "/c:" << p.component;
 }
 
 ostream & operator<< (ostream & os, sendrecv_pseudoregion_t const & srp)
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
index 033e73a05..66037bdc7 100644
--- a/Carpet/CarpetLib/src/region.hh
+++ b/Carpet/CarpetLib/src/region.hh
@@ -2,20 +2,29 @@
 #define REGION_HH
 
 #include <iostream>
+#include <vector>
 
 #include "defs.hh"
 #include "bbox.hh"
+#include "fulltree.hh"
 #include "vect.hh"
 
 
 
 // Region description
 struct region_t {
-  ibbox  extent;                            // extent
-  b2vect outer_boundaries;                  // outer boundaries
-  int    map;                               // map to which this
-                                            // region belongs
-  int    processor;                         // processor number
+  ibbox        extent;           // extent
+  b2vect       outer_boundaries; // outer boundaries
+  int          map;              // map to which this region belongs
+  int          processor;        // processor number
+  ipfulltree * processors;       // processor decomposition
+  
+  region_t ();
+  region_t (region_t const & a);
+  region_t & operator= (region_t const & a);
+  ~region_t ();
+  
+  bool invariant () const;
 };
 
 
@@ -29,54 +38,62 @@ bool operator!= (region_t const & a, region_t const & b)
 
 
 
-inline size_t memoryof (region_t const & reg)
-{
-  return
-    memoryof (reg.extent) +
-    memoryof (reg.outer_boundaries) +
-    memoryof (reg.map) +
-    memoryof (reg.processor);
-}
+void
+combine_regions (vector<region_t> const & oldregs,
+                 vector<region_t> & newregs);
+
+
+
+size_t memoryof (region_t const & reg);
 
 istream & operator>> (istream & is, region_t       & reg);
 ostream & operator<< (ostream & os, region_t const & reg);
 
 
 
-// A pseudoregion is almost a region; it is a bbox that lives on a
-// certain processor.  Pseudoregions are a compact way to store
-// information about what processors needs to send data to what other
-// processors during synchronisation or regridding.
+// A pseudoregion is almost a region; it is a bbox that belongs to a
+// certain component.  Pseudoregions are a compact way to store
+// information about what components needs to send data to what other
+// components during synchronisation or regridding.
 struct pseudoregion_t {
   ibbox extent;
-  int processor;
+  int component;
   pseudoregion_t ()
   {
   }
-  pseudoregion_t (ibbox const & extent_, int const processor_)
-    : extent (extent_), processor (processor_)
+  pseudoregion_t (ibbox const & extent_, int const component_)
+    : extent (extent_), component (component_)
   {
   }
 };
 
+bool operator== (pseudoregion_t const & a, pseudoregion_t const & b);
+inline
+bool operator!= (pseudoregion_t const & a, pseudoregion_t const & b)
+{
+  return not (a == b);
+}
+
 inline size_t memoryof (pseudoregion_t const & p)
 {
   return
     memoryof (p.extent) +
-    memoryof (p.processor);
+    memoryof (p.component);
 }
 
 ostream & operator<< (ostream & os, pseudoregion_t const & p);
 
+
+
 struct sendrecv_pseudoregion_t {
   pseudoregion_t send, recv;
   sendrecv_pseudoregion_t ()
   {
   }
-  sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_processor,
-                           ibbox const & recv_extent, int const recv_processor)
-    : send (pseudoregion_t (send_extent, send_processor)),
-      recv (pseudoregion_t (recv_extent, recv_processor))
+  sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_component,
+                           ibbox const & recv_extent,  int const recv_component)
+    : send (pseudoregion_t (send_extent, send_component)),
+      recv (pseudoregion_t (recv_extent, recv_component))
   {
   }
 };