CarpetLib: Add new datatype region_t

The new datatype region_t combines an extent (a bbox), an outer
boundary descriptor, a refinement descriptor, and a processor number:

struct region_t {
  ibbox  extent;                            // extent
  b2vect outer_boundaries;                  // outer boundaries
  b2vect refinement_boundaries;             // refinement boundaries
  int    map;                               // map to which this
                                            // region belongs
  int    processor;                         // processor number
};

These quantities are often used together, and combining them into a
single datatype simplifies the code significantly.

Adapt gh, dh, etc. to use this new datatype.

This is a major API change.

darcs-hash:20070112204130-dae7b-92cad546187b0fe499e8cfc38b2e26614a4f608c.gz

10 files changed, 210 insertions, 147 deletions

diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index 1b3da8365..98a160bf9 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: region.hh in region.hh
 includes header: vect.hh in vect.hh
 
 includes header: commstate.hh in commstate.hh
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 9b9715d47..09a078141 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -10,33 +10,13 @@
 
 #include "bbox.hh"
 #include "defs.hh"
+#include "region.hh"
 #include "vect.hh"
 
 using namespace std;
 
 
 
-istream& operator>> (istream& is, grid_structure_t& gs)
-{
-  skipws (is);
-  consume (is, '{');
-  is >> gs.bbss;
-  skipws (is);
-  consume (is, ';');
-  is >> gs.obss;
-  skipws (is);
-  consume (is, '}');
-  return is;
-}
-
-ostream& operator<< (ostream& os, grid_structure_t const& gs)
-{
-  os << "{" << gs.bbss << ";" << gs.obss << "}";
-  return os;
-}
-
-
-
 template <typename T>
 inline T ipow_helper (T x, unsigned int y)
 {
@@ -183,25 +163,30 @@ template CCTK_REAL ipow (CCTK_REAL x, int y);
 template vect<int,3> ipow (vect<int,3> x, int y);
 
 template istream& input (istream& os, vector<int>& v);
-template istream& input (istream& os, vector<grid_structure_t>& v);
 template istream& input (istream& os, vector<bbox<int,3> >& v);
 template istream& input (istream& os, vector<bbox<CCTK_REAL,3> >& v);
 template istream& input (istream& os, vector<vector<bbox<int,3> > >& v);
 template istream& input (istream& os, vector<vector<bbox<CCTK_REAL,3> > >& v);
+template istream& input (istream& os, vector<region_t>& v);
+template istream& input (istream& os, vector<vector<region_t> >& v);
+template istream& input (istream& os, vector<vector<vector<region_t> > >& v);
 template istream& input (istream& os, vector<vect<int,3> >& v);
 template istream& input (istream& os, vector<vect<vect<bool,2>,3> >& v);
 template istream& input (istream& os, vector<vector<vect<vect<bool,2>,3> > >& v);
 
 template ostream& output (ostream& os, const list<bbox<int,3> >& l);
+template ostream& output (ostream& os, const list<region_t>& l);
 template ostream& output (ostream& os, const set<bbox<int,3> >& s);
 template ostream& output (ostream& os, const set<bboxset<int,3> >& s);
 template ostream& output (ostream& os, const stack<bbox<int,3> >& 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<grid_structure_t>& v);
 template ostream& output (ostream& os, const vector<bbox<int,3> >& v);
 template ostream& output (ostream& os, const vector<bbox<CCTK_REAL,3> >& v);
+template ostream& output (ostream& os, const vector<region_t>& v);
+template ostream& output (ostream& os, const vector<vector<region_t> >& v);
+template ostream& output (ostream& os, const vector<vector<vector<region_t> > >& v);
 template ostream& output (ostream& os, const vector<list<bbox<int,3> > >& v);
 template ostream& output (ostream& os, const vector<vector<int> >& v);
 template ostream& output (ostream& os, const vector<vector<CCTK_REAL> >& v);
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 57d57a6b5..21272b58e 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -52,17 +52,6 @@ typedef vect<vect<int,dim>,2>  i2vect;
 
 
 
-// Grid structure description
-struct grid_structure_t {
-  vector <vector <ibbox> >  bbss; // refinement regions [reflevel][component]
-  vector <vector <bbvect> > obss; // outer boundaries   [reflevel][component]
-};
-
-istream& operator>> (istream& is, grid_structure_t& gs);
-ostream& operator<< (ostream& os, grid_structure_t const& gs);
-
-
-
 // A general type
 enum centering { vertex_centered, cell_centered };
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 59eff7286..8c2fe6cfd 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -112,7 +112,7 @@ void dh::allocate_bboxes ()
     for (int rl=0; rl<h.reflevels(); ++rl) {
       boxes.at(ml).at(rl).resize(h.components(rl));
       for (int c=0; c<h.components(rl); ++c) {
-        const ibbox intr = h.extents().at(ml).at(rl).at(c);
+        const ibbox intr = h.extent(ml,rl,c);
         dboxes & b = boxes.at(ml).at(rl).at(c);
 
 
@@ -130,7 +130,7 @@ void dh::allocate_bboxes ()
         i2vect odist(0);        // for owned regions
         for (int f=0; f<2; ++f) {
           for (int d=0; d<dim; ++d) {
-            if (h.outer_boundaries().at(rl).at(c)[d][f]) {
+            if (h.outer_boundaries(rl,c)[f][d]) {
               dist[f][d] = 0;
               boxes.at(ml).at(rl).at(c).is_interproc[d][f] = false;
             } else {
@@ -144,7 +144,7 @@ void dh::allocate_bboxes ()
               dist1[f][d] = is_empty ? 0 : dist[f][d];
               ibset bnd = intr.expand(dist1[0], dist1[1]) - intr;
               for (int cc=0; cc<h.components(rl); ++cc) {
-                bnd -= h.extents().at(ml).at(rl).at(cc);
+                bnd -= h.extent(ml,rl,cc);
               }
               bool const is_interproc = bnd.empty();
               boxes.at(ml).at(rl).at(c).is_interproc[d][f] = is_interproc;
@@ -170,7 +170,7 @@ void dh::allocate_bboxes ()
       
       // Make owned regions disjoint
       for (int c=0; c<h.components(rl); ++c) {
-        const ibbox intr = h.extents().at(ml).at(rl).at(c);
+        const ibbox intr = h.extent(ml,rl,c);
         dboxes & b = boxes.at(ml).at(rl).at(c);
         
         // 1. Remove all other interiors from this owned region
@@ -634,7 +634,7 @@ void dh::check_bboxes (dh::dboxes & box,
           }
           
           // Calculate boundary bboxes
-          b2vect const obs = xpose (h.outer_boundaries().at(rl).at(c));
+          b2vect const obs = h.outer_boundaries(rl,c);
           ibset bnd =
             box.interior - 
             box.interior.expand (- (ivect (obs[0]) * nboundaryzones[0]),
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index a78b35def..e81a3811f 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -80,7 +80,7 @@ void ggf::set_timelevels (const int ml, const int rl, const int new_timelevels)
       for (int tl=timelevels(ml,rl); tl<new_timelevels; ++tl) {
         storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
         storage.at(ml).at(rl).at(c).at(tl)->allocate
-          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.processor(rl,c));
       } // for tl
     } // for c
     
@@ -130,7 +130,7 @@ void ggf::recompose_allocate (const int rl)
       for (int tl=0; tl<timelevels(ml,rl); ++tl) {
         storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
         storage.at(ml).at(rl).at(c).at(tl)->allocate
-          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.processor(rl,c));
       } // for tl
     } // for c
   } // for ml
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 95b216a08..eb523551c 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -37,26 +37,18 @@ gh::gh (const vector<ivect> & reffacts_, const centering refcent_,
 gh::~gh () { }
 
 // Modifiers
-void gh::regrid (const mexts& exts,
-                 const rbnds& outer_bounds,
-                 const rprocs& procs)
+void gh::regrid (mregs const & regs)
 {
   DECLARE_CCTK_PARAMETERS;
   
   // Save the old grid hierarchy
-  _oldextents = _extents;
-  _oldouter_boundaries = _outer_boundaries;
-  _oldprocessors = _processors;
-  
-  _extents = exts;
-  _outer_boundaries = outer_bounds;
-  _processors = procs;
+  _oldregions = _regions;
+  _regions = regs;
   
   // Consistency checks
   
   // nota bene: there might be 0 refinement levels.
   
-  check_processor_number_consistency ();
   check_multigrid_consistency ();
   check_component_consistency ();
   check_base_grid_extent ();
@@ -83,8 +75,8 @@ bool gh::recompose (const int rl,
                     const bool do_prolongate)
 {
   // Handle changes in number of mglevels
-  if (_oldextents.size() != _extents.size()) {
-    _oldextents.resize (_extents.size());
+  if (_oldregions.size() != _regions.size()) {
+    _oldregions.resize (_regions.size());
   }
   
   bool const do_recompose = level_did_change(rl);
@@ -98,13 +90,10 @@ bool gh::recompose (const int rl,
     
     // Overwrite old with new grid hierarchy
     for (int ml=0; ml<mglevels(); ++ml) {
-      _oldextents.at(ml).resize (_extents.at(ml).size());
-      _oldextents.at(ml).at(rl) = _extents.at(ml).at(rl);
+      _oldregions.at(ml).resize (_regions.at(ml).size());
+      _oldregions.at(ml).at(rl) = _regions.at(ml).at(rl);
     }
-    _oldouter_boundaries.resize (_outer_boundaries.size());
-    _oldouter_boundaries.at(rl) = _outer_boundaries.at(rl);
-    _oldprocessors.resize (_processors.size());
-    _oldprocessors.at(rl) = _processors.at(rl);
+    
   }
   
   return do_recompose;
@@ -113,34 +102,21 @@ bool gh::recompose (const int rl,
 bool gh::level_did_change (const int rl) const
 {
   // Find out whether this level changed
-  if (_extents.size() != _oldextents.size()) return true;
+  if (_regions.size() != _oldregions.size()) return true;
   for (int ml=0; ml<mglevels(); ++ml) {
     assert (rl>=0 and rl<reflevels());
-    if (rl >= (int)_oldextents.at(ml).size()) return true;
-    if (_extents.at(ml).at(rl).size() != _oldextents.at(ml).at(rl).size()) {
+    if (rl >= (int)_oldregions.at(ml).size()) return true;
+    if (_regions.at(ml).at(rl).size() != _oldregions.at(ml).at(rl).size()) {
       return true;
     }
     for (int c=0; c<components(rl); ++c) {
-      if (_extents.at(ml).at(rl).at(c) != _oldextents.at(ml).at(rl).at(c)) {
+      if (_regions.at(ml).at(rl).at(c) != _oldregions.at(ml).at(rl).at(c)) {
         return true;
       }
-      if (_processors.at(rl).at(c) != _oldprocessors.at(rl).at(c)) return true;
     } // for c
   } // for ml
   return false;
 }
-
-void gh::check_processor_number_consistency ()
-{
-  for (int rl=0; rl<reflevels(); ++rl) {
-    assert (processors().size() == extents().at(0).size());
-    assert (outer_boundaries().size() == extents().at(0).size());
-    for (int c=0; c<components(rl); ++c) {
-      assert (processors().at(rl).size() == extents().at(0).at(rl).size());
-      assert (outer_boundaries().at(rl).size() == extents().at(0).at(rl).size());
-    }
-  }
-}
   
 void gh::check_multigrid_consistency ()
 {
@@ -148,14 +124,14 @@ void gh::check_multigrid_consistency ()
   for (int ml=1; ml<mglevels(); ++ml) {
     for (int rl=0; rl<reflevels(); ++rl) {
       for (int c=0; c<components(rl); ++c) {
-	assert (all(extents().at(ml).at(rl).at(c).stride()
-		    == ivect(mgfact) * extents().at(ml-1).at(rl).at(c).stride()));
+	assert (all(extent(ml,rl,c).stride()
+		    == ivect(mgfact) * extent(ml-1,rl,c).stride()));
         // TODO: put the check back in, taking outer boundaries into
         // account
 #if 0
-	assert (extents().at(ml).at(rl).at(c)
-		.contracted_for(extents().at(ml-1).at(rl).at(c))
-		.is_contained_in(extents().at(ml-1).at(rl).at(c)));
+	assert (extent(ml,rl,c)
+		.contracted_for(extent(ml-1,rl,c))
+		.is_contained_in(extent(ml-1,rl,c)));
 #endif
       }
     }
@@ -168,12 +144,12 @@ void gh::check_component_consistency ()
     for (int rl=0; rl<reflevels(); ++rl) {
       assert (components(rl)>0);
       for (int c=0; c<components(rl); ++c) {
-        const ibbox &b  = extents().at(ml).at(rl).at(c);
-        const ibbox &b0 = extents().at(ml).at(rl).at(0);
+        const ibbox &b  = extent(ml,rl,c);
+        const ibbox &b0 = extent(ml,rl,0);
 	assert (all(b.stride() == b0.stride()));
 	assert (b.is_aligned_with(b0));
         for (int cc=c+1; cc<components(rl); ++cc) {
-          assert ((b & extents().at(ml).at(rl).at(cc)).empty());
+          assert ((b & extent(ml,rl,cc)).empty());
         }
       }
     }
@@ -187,7 +163,7 @@ void gh::check_base_grid_extent ()
       // TODO: put the check back in, taking outer boundaries into
       // account
 #if 0
-      assert (extents().at(0).at(c).at(0).is_contained_in(baseextent));
+      assert (extent(0,c,0).is_contained_in(baseextent));
 #endif
     }
   }
@@ -200,29 +176,28 @@ void gh::check_refinement_levels ()
   bool have_error = false;
   for (int ml=0; ml<mglevels(); ++ml) {
     for (int rl=1; rl<reflevels(); ++rl) {
-      assert (all(extents().at(ml).at(rl-1).at(0).stride()
+      assert (all(extent(ml,rl-1,0).stride()
                   == ((reffacts.at(rl) / reffacts.at(rl-1))
-                      * extents().at(ml).at(rl).at(0).stride())));
+                      * extent(ml,rl,0).stride())));
       // Check contained-ness:
       // first take all coarse grids ...
       ibset all;
       for (int c=0; c<components(rl-1); ++c) {
-        all |= extents().at(ml).at(rl-1).at(c);
+        all |= extent(ml,rl-1,c);
       }
 #if 0
       // ... remember their size ...
       const int sz = all.size();
       // ... then add the coarsified fine grids ...
       for (int c=0; c<components(rl); ++c) {
-        all |= extents().at(ml).at(rl).at(c).contracted_for(extents().at(ml).at(rl-1).at(0));
+        all |= extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
       }
       // ... and then check the sizes:
       assert (all.size() == sz);
 #endif
       for (int c=0; c<components(rl); ++c) {
         ibset const finebox
-          = (extents().at(ml).at(rl).at(c).contracted_for
-             (extents().at(ml).at(rl-1).at(0)));
+          = (extent(ml,rl,c).contracted_for (extent(ml,rl-1,0)));
         if (! (finebox <= all)) {
           if (! have_error) {
             cout << "The following components are not properly nested, i.e., they are not" << endl
@@ -240,7 +215,7 @@ void gh::check_refinement_levels ()
       for (int rl=0; rl<reflevels(); ++rl) {
         for (int c=0; c<components(rl); ++c) {
           cout << "   ml " << ml << " rl " << rl << " c " << c << ":   "
-               << extents().at(ml).at(rl).at(c) << endl;
+               << extent(ml,rl,c) << endl;
         }
       }
     }
@@ -257,7 +232,7 @@ void gh::calculate_base_extents_of_all_levels ()
       _bases.at(ml).at(rl) = ibbox();
       ibbox &bb = _bases.at(ml).at(rl);
       for (int c=0; c<components(rl); ++c) {
-        bb = bb.expanded_containing(extents().at(ml).at(rl).at(c));
+        bb = bb.expanded_containing(extent(ml,rl,c));
       }
     }
   }
@@ -308,9 +283,10 @@ void gh::do_output_bboxes (ostream& os) const
         os << endl;
         os << "gh bboxes:" << endl;
         os << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
-        os << "extent=" << extents().at(ml).at(rl).at(c) << endl;
-        os << "outer_boundary=" << outer_boundaries().at(rl).at(c) << endl;
-        os << "processor=" << processors().at(rl).at(c) << endl;
+        os << "extent=" << extent(ml,rl,c) << endl;
+        os << "outer_boundaries=" << outer_boundaries(rl,c) << endl;
+        os << "refinement_boundaries=" << refinement_boundaries(rl,c) << endl;
+        os << "processor=" << processor(rl,c) << endl;
       }
     }
   }
@@ -335,9 +311,7 @@ ostream& gh::output (ostream& os) const
   os << "gh:"
      << "reffacts=" << reffacts << ",refcentering=" << refcent << ","
      << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
-     << "extents=" << extents() << ","
-     << "outer_boundaries=" << outer_boundaries() << ","
-     << "processors=" << processors() << ","
+     << "regions=" << regions() << ","
      << "dhs={";
   const char * sep = "";
   for (list<dh*>::const_iterator d = dhs.begin();
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index a3e943785..cc9d148f2 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -10,6 +10,7 @@
 #include "bboxset.hh"
 #include "defs.hh"
 #include "dist.hh"
+#include "region.hh"
 #include "vect.hh"
 
 using namespace std;
@@ -33,15 +34,9 @@ class gh {
 public:
   
   // Types
-  typedef vector<ibbox> cexts;	// ... for each component
-  typedef vector<cexts> rexts;	// ... for each refinement level
-  typedef vector<rexts> mexts;	// ... for each multigrid level
-  
-  typedef vector<bbvect> cbnds;	// ... for each component
-  typedef vector<cbnds> rbnds;	// ... for each refinement level
-  
-  typedef vector<int>    cprocs; // ... for each component
-  typedef vector<cprocs> rprocs; // ... for each refinement level
+  typedef vector<region_t> cregs; // ... for each component
+  typedef vector<cregs> rregs;    // ... for each refinement level
+  typedef vector<rregs> mregs;    // ... for each multigrid level
   
 public:				// should be readonly
   
@@ -54,18 +49,14 @@ public:				// should be readonly
   
   const ibbox baseextent;
   
-  
 private:
   vector<vector<ibbox> > _bases; // [ml][rl]
-  // TODO: invent structure for this
-  mexts _extents;		// extents of all grids
-  rbnds _outer_boundaries;	// boundary descriptions of all grids
-  rprocs _processors;		// processor numbers of all grids
-
-  mexts _oldextents;            // a copy, used during regridding
-  rbnds _oldouter_boundaries;
-  rprocs _oldprocessors;
-
+  
+  // Extents and properties of all grids
+  mregs _regions;
+  // A copy, used during regridding
+  mregs _oldregions;
+  
   list<th*> ths;		// list of all time hierarchies
   list<dh*> dhs;		// list of all data hierarchies
   
@@ -80,65 +71,64 @@ public:
   virtual ~gh ();
   
   // Modifiers
-  void regrid (const mexts& exts,
-               const rbnds& outer_bounds,
-               const rprocs& procs);
+  void regrid (mregs const & regs);
   bool recompose (const int rl,
                   const bool do_prolongate);
 private:
   bool level_did_change (const int rl) const;
   
+  // Accessors
+  
 public:
-  const mexts & extents() const
+  mregs const & regions () const
   {
-    return _extents;
+    return _regions;
   }
 
-  const rbnds & outer_boundaries() const
+  const vector<vector<ibbox> > & bases() const
+  {
+    return _bases;
+  }
+  
+  ibbox extent (const int m, const int rl, const int c) const
   {
-    return _outer_boundaries;
+    return _regions.at(m).at(rl).at(c).extent;
+  }
+  
+  b2vect outer_boundaries (const int rl, const int c) const
+  {
+    return _regions.at(0).at(rl).at(c).outer_boundaries;
   }
 
-  const rprocs & processors() const
+  b2vect refinement_boundaries (const int rl, const int c) const
   {
-    return _processors;
+    return _regions.at(0).at(rl).at(c).refinement_boundaries;
   }
 
-  const vector<vector<ibbox> > & bases() const
+  int processor (const int rl, const int c) const
   {
-    return _bases;
+    return _regions.at(0).at(rl).at(c).processor;
   }
-  
-  // Accessors
+
   int mglevels () const
   {
-    return (int)_extents.size();
+    return (int)_regions.size();
   }
   
   int reflevels () const
   {
     if (mglevels() == 0) return 0;
-    return (int)_extents.at(0).size();
+    return (int)_regions.at(0).size();
   }
   
   int components (const int rl) const
   {
-    return (int)_extents.at(0).at(rl).size();
-  }
-  
-  bbvect outer_boundary (const int rl, const int c) const
-  {
-    return _outer_boundaries.at(rl).at(c);
-  }
-  
-  int proc (const int rl, const int c) const
-  {
-    return _processors.at(rl).at(c);
+    return (int)_regions.at(0).at(rl).size();
   }
 
   bool is_local (const int rl, const int c) const
   {
-    return proc(rl,c) == dist::rank();
+    return processor(rl,c) == dist::rank();
   }
   
   int local_components (const int rl) const;
@@ -155,7 +145,6 @@ public:
   virtual ostream& output (ostream& os) const;
 
 private:
-  void check_processor_number_consistency ();
   void check_multigrid_consistency ();
   void check_component_consistency ();
   void check_base_grid_extent ();
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 7a846f845..6c85dc167 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -13,6 +13,7 @@ SRCS =	bbox.cc					\
 	ggf.cc					\
 	gh.cc					\
 	mem.cc					\
+	region.cc				\
 	th.cc					\
 	vect.cc					\
 	checkindex.c				\
diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc
new file mode 100644
index 000000000..db9ef302e
--- /dev/null
+++ b/Carpet/CarpetLib/src/region.cc
@@ -0,0 +1,86 @@
+#include <iostream>
+
+#include "defs.hh"
+#include "region.hh"
+
+using namespace std;
+
+
+
+bool
+operator== (region_t const & a, region_t const & b)
+{
+  return
+    a.extent == b.extent and
+    all (all (a.outer_boundaries == b.outer_boundaries)) and
+    all (all (a.refinement_boundaries == b.refinement_boundaries)) and
+    a.map == b.map and
+    a.processor == b.processor;
+}
+
+
+
+istream &
+operator>> (istream & is, region_t & reg)
+{
+  skipws (is);
+  consume (is, "region_t");
+  skipws (is);
+  consume (is, '(');
+  
+  skipws (is);
+  consume (is, "extent");
+  skipws (is);
+  consume (is, '=');
+  is >> reg.extent;
+  skipws (is);
+  consume (is, ',');
+  
+  skipws (is);
+  consume (is, "outer_boundaries");
+  skipws (is);
+  consume (is, '=');
+  is >> reg.outer_boundaries;
+  skipws (is);
+  consume (is, ',');
+  
+  skipws (is);
+  consume (is, "refinement_boundaries");
+  skipws (is);
+  consume (is, '=');
+  is >> reg.refinement_boundaries;
+  skipws (is);
+  consume (is, ',');
+  
+  skipws (is);
+  consume (is, "map");
+  skipws (is);
+  consume (is, '=');
+  is >> reg.map;
+  skipws (is);
+  consume (is, ',');
+  
+  skipws (is);
+  consume (is, "processor");
+  skipws (is);
+  consume (is, '=');
+  is >> reg.processor;
+  skipws (is);
+  consume (is, ')');
+  
+  return is;
+}
+
+
+
+ostream &
+operator<< (ostream & os, region_t const & reg)
+{
+  os << "region_t("
+     << "extent=" << reg.extent << ","
+     << "outer_boundaries=" << reg.outer_boundaries << ","
+     << "refinement_boundaries=" << reg.refinement_boundaries << ","
+     << "map=" << reg.map << ","
+     << "processor=" << reg.processor << ")";
+  return os;
+}
diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh
new file mode 100644
index 000000000..73bc48ea0
--- /dev/null
+++ b/Carpet/CarpetLib/src/region.hh
@@ -0,0 +1,38 @@
+#ifndef REGION_HH
+#define REGION_HH
+
+#include <iostream>
+
+#include "defs.hh"
+#include "bbox.hh"
+#include "vect.hh"
+
+
+
+// Region description
+struct region_t {
+  ibbox  extent;                            // extent
+  b2vect outer_boundaries;                  // outer boundaries
+  b2vect refinement_boundaries;             // refinement boundaries
+  int    map;                               // map to which this
+                                            // region belongs
+  int    processor;                         // processor number
+};
+
+
+
+bool operator== (region_t const & a, region_t const & b);
+inline
+bool operator!= (region_t const & a, region_t const & b)
+{
+  return not operator== (a, b);
+}
+
+
+
+istream & operator>> (istream & is, region_t       & reg);
+ostream & operator<< (ostream & os, region_t const & reg);
+
+
+
+#endif // #ifndef REGION_HH