CarpetLib: Add new bboxset functions

New functions:
- construct a bboxset from a container of bboxes
- expand or contract bboxsets

5 files changed, 177 insertions, 95 deletions

diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 71fb0c9a5..7ffd61699 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -49,7 +49,27 @@ bboxset<T,D>::bboxset (const vector<list<box> >& vlb) {
   }
 }
 
-template<class T, int D>
+template<typename T, int D>
+template<typename U>
+bboxset<T,D>::bboxset (const vector<U>& vb, const bbox<T,D> U::* const v) {
+  for (typename vector<U>::const_iterator
+         vi = vb.begin(); vi != vb.end(); ++ 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) {
+  for (typename vector<U>::const_iterator
+         vi = vb.begin(); vi != vb.end(); ++ vi)
+  {
+    *this |= (*vi).*v;
+  }
+}
+
+template<typename T, int D>
 bboxset<T,D> bboxset<T,D>::poison () {
   return bboxset (bbox<T,D>::poison());
 }
@@ -427,6 +447,53 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s) {
 
 
 
+template<typename T, int D>
+typename bboxset<T,D>::box bboxset<T,D>::container () const {
+  box b;
+  for (const_iterator bi=begin(); bi!=end(); ++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;
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const {
+  // We don't know (yet?) how to shrink a set
+  assert (all (lo>=0 and hi>=0));
+  bboxset res;
+  for (const_iterator bi=begin(); bi!=end(); ++bi) {
+    res |= (*bi).expand(lo,hi);
+  }
+  return res;
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::expanded_for (const box& b) const {
+  bboxset res;
+  for (const_iterator bi=begin(); bi!=end(); ++bi) {
+    res |= (*bi).expanded_for(b);
+  }
+  return res;
+}
+
+template<typename T, int D>
+bboxset<T,D> bboxset<T,D>::contracted_for (const box& b) const {
+  bboxset res;
+  for (const_iterator bi=begin(); bi!=end(); ++bi) {
+    res |= (*bi).contracted_for(b);
+  }
+  return res;
+}
+
+
+
 // Equality
 template<typename T, int D>
 bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const {
@@ -512,3 +579,13 @@ ostream& bboxset<T,D>::output (ostream& os) const {
 
 
 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 85eb2b2c8..5624ef4b1 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.hh
@@ -63,6 +63,10 @@ public:
   
   bboxset (const list<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 ();
   
@@ -133,6 +137,23 @@ public:
   // friend bboxset operator- <T,D>(const box& b, const bboxset& s);
   static bboxset minus (const box& b, const bboxset& s);
   
+  /** Find a bbox containing the whole set.  */
+  box container () const;
+  /** Find the pseudo-inverse.  */
+  bboxset pseudo_inverse (const int n) const;
+  
+  /** Expand (enlarge) the bbox by multiples of the stride.  */
+  bboxset expand (const vect<T,D>& lo, const vect<T,D>& hi) const;
+  bboxset expand (const vect<vect<T,D>,2>& lohi) const
+  { return expand (lohi[0], lohi[1]); }
+  
+  /** Find the smallest b-compatible box around this bbox.
+      ("compatible" means having the same stride.)  */
+  bboxset expanded_for (const box& b) const;
+  
+  /** Find the largest b-compatible box inside this bbox.  */
+  bboxset contracted_for (const box& b) const;
+  
   // Equality
   bool operator== (const bboxset& s) const;
   bool operator!= (const bboxset& s) const;
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index a2b366ea2..1ba4f302d 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -330,6 +330,8 @@ template CCTK_REAL ipow (CCTK_REAL 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<dh*> const & l);
@@ -342,6 +344,7 @@ 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<ivect> const & v);
 template size_t memoryof (vector<i2vect> const & v);
 template size_t memoryof (vector<fulltree <int,dim,pseudoregion_t> *> const & f);
@@ -367,6 +370,7 @@ 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<ivect>& v);
 template istream& input (istream& os, vector<bbvect>& v);
 template istream& input (istream& os, vector<i2vect>& v);
@@ -394,6 +398,7 @@ 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<rbbox>& v);
+template ostream& output (ostream& os, const vector<ibset>& v);
 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);
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 51ccc63ef..b1f4e3937 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -444,12 +444,7 @@ regrid (bool const do_init)
       ibbox const domain_enlarged = domain_active.expand (safedist);
       
       // All owned regions
-      ibset allowned;
-      for (int c = 0; c < h.components(rl); ++ c) {
-        full_dboxes const & box = full_level.AT(c);
-        allowned += box.owned;
-      }
-      allowned.normalize();
+      ibset const allowned (full_level, & full_dboxes::owned);
       ASSERT_rl (allowned <= domain_active,
                  "The owned regions must be contained in the active part of the domain");
       
@@ -492,14 +487,7 @@ regrid (bool const do_init)
       
       
       // The conjunction of all buffer zones must equal allbuffers
-      // cerr << "QQQ: regrid[d]" << endl;
-      
-      ibset allbuffers1;
-      for (int c = 0; c < h.components(rl); ++ c) {
-        full_dboxes const & box = full_level.AT(c);
-        allbuffers1 += box.buffers;
-      }
-      allbuffers1.normalize();
+      ibset const allbuffers1 (full_level, & full_dboxes::buffers);
       ASSERT_rl (allbuffers1 == allbuffers,
                  "Buffer zone consistency check");
       
@@ -574,17 +562,9 @@ regrid (bool const do_init)
           
           ibset needrecv = box.active;
           
-          ibset contracted_oactive;
-          for (ibset::const_iterator
-                 ai = obox.active.begin(); ai != obox.active.end(); ++ ai)
-          {
-            ibbox const & oactive = * ai;
-            contracted_oactive += oactive.contracted_for (box.interior);
-          }
-          contracted_oactive.normalize();
-          
-          ibset ovlp = needrecv & contracted_oactive;
-          ovlp.normalize();
+          ibset const contracted_oactive
+            (obox.active.contracted_for (box.interior));
+          ibset const ovlp = needrecv & contracted_oactive;
           
           for (ibset::const_iterator
                  ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
@@ -628,17 +608,8 @@ regrid (bool const do_init)
           
           i2vect const stencil_size = i2vect (prolongation_stencil_size(rl));
           
-          ibset expanded_active;
-          for (ibset::const_iterator
-                 ai = box.active.begin(); ai != box.active.end(); ++ ai)
-          {
-            ibbox const & active = * ai;
-            expanded_active += active.expanded_for (obox.interior);
-          }
-          expanded_active.normalize();
-          
-          ibset ovlp = oneedrecv & expanded_active;
-          ovlp.normalize();
+          ibset const expanded_active (box.active.expanded_for (obox.interior));
+          ibset const ovlp = oneedrecv & expanded_active;
           
           for (ibset::const_iterator
                  ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
@@ -688,19 +659,15 @@ regrid (bool const do_init)
           for (int cc = 0; cc < h.components(orl); ++ cc) {
             full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc);
             
-            ibset contracted_oactive;
-            for (ibset::const_iterator
-                   ai = obox.active.begin(); ai != obox.active.end(); ++ ai)
-            {
-              ibbox const & oactive = * ai;
-              // untested for cell centering
-              contracted_oactive +=
-                oactive.contracted_for (box.interior).expand (reffact);
-            }
-            contracted_oactive.normalize();
-            
-            ibset ovlp = needrecv & contracted_oactive;
-            ovlp.normalize();
+#if 0
+            // untested for cell centering
+            ibset const expanded_oactive
+              (obox.active.contracted_for (box.interior).expand (reffact));
+#else
+            ibset const expanded_oactive
+              (obox.active.expanded_for (box.interior).expand (reffact));
+#endif
+            ibset const ovlp = needrecv & expanded_oactive;
             
             for (ibset::const_iterator
                    ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
@@ -846,19 +813,15 @@ regrid (bool const do_init)
           for (int cc = 0; cc < h.components(orl); ++ cc) {
             full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc);
             
-            ibset contracted_oactive;
-            for (ibset::const_iterator
-                   ai = obox.active.begin(); ai != obox.active.end(); ++ ai)
-            {
-              ibbox const & oactive = * ai;
-              // untested for cell centering
-              contracted_oactive +=
-                oactive.contracted_for (box.interior).expand (reffact);
-            }
-            contracted_oactive.normalize();
-            
-            ibset ovlp = needrecv & contracted_oactive;
-            ovlp.normalize();
+#if 0
+            // untested for cell centering
+            ibset const expanded_oactive
+              (obox.active.contracted_for (box.interior).expand (reffact));
+#else
+            ibset const expanded_oactive
+              (obox.active.expanded_for (box.interior).expand (reffact));
+#endif
+            ibset const ovlp = needrecv & expanded_oactive;
             
             for (ibset::const_iterator
                    ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
@@ -917,29 +880,14 @@ regrid (bool const do_init)
             // Refinement restriction may fill all active points, and
             // must use all active points
             
-            ibset needrecv;
-            for (ibset::const_iterator
-                   ai = box.active.begin(); ai != box.active.end(); ++ ai)
-            {
-              ibbox const & active = * ai;
-              needrecv += active.contracted_for (obox0.interior);
-            }
-            needrecv.normalize();
+            ibset needrecv (box.active.contracted_for (obox0.interior));
             
             for (int cc = 0; cc < h.components(orl); ++ cc) {
               full_dboxes & obox = full_boxes.AT(ml).AT(orl).AT(cc);
               
-              ibset contracted_active;
-              for (ibset::const_iterator
-                     ai = box.active.begin(); ai != box.active.end(); ++ ai)
-              {
-                ibbox const & active = * ai;
-                contracted_active += active.contracted_for (obox0.interior);
-              }
-              contracted_active.normalize();
-              
-              ibset ovlp = obox.active & contracted_active;
-              ovlp.normalize();
+              ibset const contracted_active
+                (box.active.contracted_for (obox0.interior));
+              ibset const ovlp = obox.active & contracted_active;
               
               for (ibset::const_iterator
                      ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
@@ -1066,19 +1014,15 @@ regrid (bool const do_init)
             for (int cc = 0; cc < h.components(orl); ++ cc) {
               full_dboxes const & obox = full_boxes.AT(ml).AT(orl).AT(cc);
               
-              ibset contracted_oactive;
-              for (ibset::const_iterator
-                     ai = obox.active.begin(); ai != obox.active.end(); ++ ai)
-              {
-                ibbox const & oactive = * ai;
-                // untested for cell centering
-                contracted_oactive +=
-                  oactive.contracted_for (box.interior).expand (reffact);
-              }
-              contracted_oactive.normalize();
-              
-              ibset ovlp = needrecv & contracted_oactive;
-              ovlp.normalize();
+#if 0
+              // untested for cell centering
+              ibset const expanded_oactive
+                (obox.active.contracted_for (box.interior).expand (reffact));
+#else
+              ibset const expanded_oactive
+                (obox.active.expanded_for (box.interior).expand (reffact));
+#endif
+              ibset const ovlp = needrecv & expanded_oactive;
               
               for (ibset::const_iterator
                      ri = ovlp.begin(); ri != ovlp.end(); ++ ri)
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc
index 4a573d6af..b02a46979 100644
--- a/Carpet/CarpetLib/src/vect.cc
+++ b/Carpet/CarpetLib/src/vect.cc
@@ -4,6 +4,7 @@
 #include "cctk.h"
 
 #include "defs.hh"
+#include "bboxset.hh"
 
 #include "vect.hh"
 
@@ -83,3 +84,37 @@ template void vect<vect<bool,2>,dim>::output (ostream& os) const;
 template void vect<vect<int,2>,dim>::output (ostream& os) const;
 template void vect<vect<bool,dim>,2>::output (ostream& os) const;
 template void vect<vect<int,dim>,2>::output (ostream& os) const;
+template void vect<vect<CCTK_REAL,dim>,2>::output (ostream& os) const;
+
+
+
+// Instantiate for bboxset class
+
+#define DEFINE_FAKE_VECT_OPERATIONS(T,D)                                \
+template<> vect<T,D> vect<T,D>::dir (const int d) { assert(0); }        \
+template<> vect<T,D> vect<T,D>::seq () { assert(0); }                   \
+template<> vect<T,D> vect<T,D>::seq (const int n) { assert(0); }        \
+template<> vect<T,D> vect<T,D>::seq (const int n, const int s) { assert(0); } \
+template<> vect<T,D>& vect<T,D>::operator*= (const vect<T,D>&) { assert(0); } \
+template<> vect<T,D>& vect<T,D>::operator*= (const T&) { assert(0); }   \
+template<> vect<T,D>& vect<T,D>::operator/= (const vect<T,D>&) { assert(0); } \
+template<> vect<T,D>& vect<T,D>::operator/= (const T&) { assert(0); }   \
+template<> vect<T,D>& vect<T,D>::operator%= (const vect<T,D>&) { assert(0); } \
+template<> vect<T,D>& vect<T,D>::operator%= (const T&) { assert(0); }   \
+template<> vect<T,D>& vect<T,D>::operator^= (const vect<T,D>&) { assert(0); } \
+template<> vect<T,D>& vect<T,D>::operator^= (const T&) { assert(0); }   \
+template<> vect<T,D> vect<T,D>::operator+ () const { assert(0); }       \
+template<> vect<T,D> vect<T,D>::operator- () const { assert(0); }       \
+template<> vect<T,D> vect<T,D>::operator~ () const { assert(0); }       \
+template class vect<T,D>;                                               \
+template size_t memoryof (const vect<T,D>&);                            \
+template istream& operator>> (istream& is, vect<T,D>&);                 \
+template ostream& operator<< (ostream& os, const vect<T,D>&);
+
+typedef bboxset<int,dim> T1;
+typedef vect<bboxset<int,dim>,2> T2;
+
+DEFINE_FAKE_VECT_OPERATIONS(T1,dim)
+DEFINE_FAKE_VECT_OPERATIONS(T2,dim)
+
+#undef DEFINE_FAKE_VECT_OPERATIONS