diff options
Diffstat (limited to 'Carpet/CarpetLib/src')
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.cc | 793 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.hh | 383 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset1.cc | 55 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset1.hh | 1189 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset2.cc | 38 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset2.hh | 1369 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/defs.cc | 144 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/defs.hh | 38 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/make.code.defn | 2 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/region.hh | 1 |
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" |