From 9791222a06c996805b3509be67f30aa731fcae10 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 26 May 2013 16:25:12 -0400 Subject: CarpetLib: New class bboxset2 Rename bboxset to bboxset1. Implement new class bboxset2, which uses a different internal datastructure than bboxset1. Make bboxset a typedef for either bboxset1 (default) or bboxset2, as selected by the compile-time macro CARPET_BBOXSET2. Disable all bboxset2 code if CARPET_NO_BBOXSET2 is given, since bboxset2 uses newer C++ constructs not available on some older compilers. --- Carpet/CarpetLib/src/bboxset.cc | 793 -------------------- Carpet/CarpetLib/src/bboxset.hh | 383 +--------- Carpet/CarpetLib/src/bboxset1.cc | 55 ++ Carpet/CarpetLib/src/bboxset1.hh | 1189 ++++++++++++++++++++++++++++++ Carpet/CarpetLib/src/bboxset2.cc | 38 + Carpet/CarpetLib/src/bboxset2.hh | 1369 +++++++++++++++++++++++++++++++++++ Carpet/CarpetLib/src/defs.cc | 144 ++-- Carpet/CarpetLib/src/defs.hh | 38 +- Carpet/CarpetLib/src/make.code.defn | 2 + Carpet/CarpetLib/src/region.hh | 1 + 10 files changed, 2786 insertions(+), 1226 deletions(-) create mode 100644 Carpet/CarpetLib/src/bboxset1.cc create mode 100644 Carpet/CarpetLib/src/bboxset1.hh create mode 100644 Carpet/CarpetLib/src/bboxset2.cc create mode 100644 Carpet/CarpetLib/src/bboxset2.hh (limited to 'Carpet/CarpetLib') 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 -#include -#include -#include -#include -#include - -#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 -bboxset::bboxset () - : skip_normalize(false) -{ - assert (invariant()); -} - -template -bboxset::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 -bboxset::bboxset (const bboxset& s) - : bs(s.bs), skip_normalize(false) -{ - assert (invariant()); -} - -template -bboxset::bboxset (const list& lb) - : skip_normalize(false) -{ - SKIP_NORMALIZE(*this); - for (typename list::const_iterator - li = lb.begin(), le = lb.end(); li != le; ++ li) - { - *this |= *li; - } -} - -template -bboxset::bboxset (const vector& lb) - : skip_normalize(false) -{ - SKIP_NORMALIZE(*this); - for (typename vector::const_iterator - vi = lb.begin(), ve = lb.end(); vi != ve; ++ vi) - { - *this |= *vi; - } -} - -template -bboxset::bboxset (const vector >& vlb) - : skip_normalize(false) -{ - SKIP_NORMALIZE(*this); - for (typename vector >::const_iterator - vli = vlb.begin(), vle = vlb.end(); vli != vle; ++ vli) - { - *this |= bboxset (*vli); - } -} - -template -template -bboxset::bboxset (const vector& vb, const bbox U::* const v) - : skip_normalize(false) -{ - SKIP_NORMALIZE(*this); - for (typename vector::const_iterator - vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi) - { - *this |= (*vi).*v; - } -} - -template -template -bboxset::bboxset (const vector& vb, const bboxset U::* const v) - : skip_normalize(false) -{ - SKIP_NORMALIZE(*this); - for (typename vector::const_iterator - vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi) - { - *this |= (*vi).*v; - } -} - -template -bboxset bboxset::poison () { - return bboxset (bbox::poison()); -} - - - -// Invariant -template -bool bboxset::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 -void bboxset::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 buf; - typedef vector 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; dsize(); - - // Can't use operators on *this since these would call normalize again - // assert (*this == oldbs); - assert (newsize == oldsize); -} - - - -// Accessors -// cost: O(n) -template -typename bboxset::size_type bboxset::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::max() - bsz >= s); - s += bsz; - } - return s; -} - - - -// Queries - -// Containment -// cost: O(n) -template -bool bboxset::contains (const vect& 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 -bool bboxset::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 -bboxset& bboxset::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 -bboxset& bboxset::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 -bboxset& bboxset::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 -bboxset bboxset::operator+ (const box& b) const -{ - bboxset r(*this); - r += b; - assert (r.invariant()); - return r; -} - -// cost: O(n) -template -bboxset bboxset::operator+ (const bboxset& s) const -{ - bboxset r(*this); - r += s; - assert (r.invariant()); - return r; -} - - - -// Union -// cost: O(n) -template -bboxset& bboxset::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 -bboxset& bboxset::operator|= (const bboxset& s) -{ - bboxset tmp = s - *this; - add_transfer (tmp); - assert (invariant()); - return *this; -} - -// cost: O(n) -template -bboxset bboxset::operator| (const box& b) const -{ - bboxset r(*this); - r |= b; - assert (r.invariant()); - return r; -} - -// cost: O(n^2) -template -bboxset bboxset::operator| (const bboxset& s) const -{ - bboxset r(*this); - r |= s; - assert (r.invariant()); - return r; -} - - - -// Intersection -// cost: O(n) -template -bboxset bboxset::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 -bboxset bboxset::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 -bboxset& bboxset::operator&= (const box& b) -{ - *this = *this & b; - assert (invariant()); - return *this; -} - -// cost: O(n) -template -bboxset& bboxset::operator&= (const bboxset& s) -{ - *this = *this & s; - assert (invariant()); - return *this; -} - - - -// Difference -// cost: O(1) -template -bboxset bboxset::minus (const bbox& b1, const bbox& b2) -{ - assert (b1.is_aligned_with(b2)); - if (b1.empty()) return bboxset(); - if (not b1.intersects(b2)) return bboxset(b1); - vect const b2lo = max (b1.lower(), b2.lower()); - vect const b2up = min (b1.upper(), b2.upper()); - vect const & b1lo = b1.lower(); - vect const & b1up = b1.upper(); - vect const & str = b1.stride(); - bboxset r; - SKIP_NORMALIZE(r); - { - for (int d=0; d lb, ub; - bbox b; - for (int dd=0; ddd) { - lb[dd] = b1lo[dd]; - ub[dd] = b1up[dd]; - } - } - lb[d] = b1lo[d]; - ub[d] = b2lo[d] - str[d]; - b = bbox(lb,ub,str); - r += b; - lb[d] = b2up[d] + str[d]; - ub[d] = b1up[d]; - b = bbox(lb,ub,str); - r += b; - } - } - assert (r.invariant()); - return r; -} - -// cost: O(n) -template -bboxset bboxset::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 -bboxset& bboxset::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 -bboxset bboxset::operator- (const bboxset& s) const -{ - bboxset r(*this); - r -= s; - assert (r.invariant()); - return r; -} - -// cost: O(n^2) -template -bboxset bboxset::minus (const bbox& b, const bboxset& s) -{ - bboxset r = bboxset(b) - s; - assert (r.invariant()); - return r; -} - - - -// cost: O(n) -template -typename bboxset::box bboxset::container () const -{ - box b; - for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) { - b = b.expanded_containing(*bi); - } - return b; -} - -template -bboxset bboxset::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 -bboxset bboxset::expand (const vect& lo, const vect& 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 -bboxset bboxset::expand (const vect& lo, const vect& hi, - const vect& denom) const -{ - assert (all(denom > vect(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 -bboxset bboxset::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 -bboxset bboxset::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 -bool bboxset::operator<= (const bboxset& s) const -{ - return (*this - s).empty(); -} - -template -bool bboxset::operator< (const bboxset& s) const -{ - return (*this - s).empty() && not (s - *this).empty(); -} - -template -bool bboxset::operator>= (const bboxset& s) const -{ - return s <= *this; -} - -template -bool bboxset::operator> (const bboxset& s) const -{ - return s < *this; -} - -template -bool bboxset::operator== (const bboxset& s) const -{ - return (*this <= s) && (*this >= s); -} - -template -bool bboxset::operator!= (const bboxset& s) const -{ - return not (*this == s); -} - - - -// Input -template -istream& bboxset::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 -ostream& bboxset::output (ostream& os) const -{ - T Tdummy; - os << "bboxset<" << typestring(Tdummy) << "," << D << ">:{" - << "size=" << size() << "," - << "setsize=" << setsize() << "," - << "set=" << bs - << "}"; - return os; -} - - - -template class bboxset; -template size_t memoryof (const bboxset& s); -template istream& operator>> (istream& is, bboxset& s); -template ostream& operator<< (ostream& os, const bboxset& s); - -#include "dh.hh" -#include "region.hh" - -template bboxset::bboxset (const vector& vb, const bbox dh::full_dboxes::* const v); -template bboxset::bboxset (const vector& vb, const bboxset dh::full_dboxes::* const v); -template bboxset::bboxset (const vector& vb, const bbox 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 -#include -#include -#include -#include - -#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 class bboxset; - -// template -// bboxset operator+ (const bbox& b1, const bbox& b2); -// template -// bboxset operator+ (const bbox& b, const bboxset& s); - -// template -// bboxset operator- (const bbox& b1, const bbox& b2); -// template -// bboxset operator- (const bbox& b, const bboxset& s); - -// Input -template -istream& operator>> (istream& is, bboxset& s); - -// Output -template -ostream& operator<< (ostream& os, const bboxset& s); - - - -// Bounding box set class -template -class bboxset { - - // Cost annotations depend on the number of bset elements n, and - // assume that normalization is skipped. - - struct skip_normalize_t { - bboxset& s; - bool const saved_skip_normalize; - skip_normalize_t(bboxset& 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 box; -#ifdef BBOXSET_SET - //S typedef set bset; -#endif -#ifdef BBOXSET_LIST - typedef list bset; -#endif -#ifdef BBOXSET_VECTOR - typedef vector 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& lb); - bboxset (const vector& lb); - bboxset (const vector >& vlb); - template - bboxset (const vector& vb, const bbox U::* const v); - template - bboxset (const vector& 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& 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- (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- (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& lo, const vect& hi) - const CCTK_MEMBER_ATTRIBUTE_PURE; - bboxset expand (const vect,2>& lohi) const - { return expand (lohi[0], lohi[1]); } - - /** Shift the bboxset by multiples of the stride. */ - bboxset shift (const vect& 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& lo, const vect& hi, - const vect& denom) const CCTK_MEMBER_ATTRIBUTE_PURE; - // cost: O(n^2) in general, but only O(n) for shifting - bboxset expand (const vect,2>& lohi, const vect& 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& v, const vect& 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 -inline bboxset operator+ (const bbox& b1, const bbox& b2) { - return bboxset(b1) + bboxset(b2); -} - -template -inline bboxset operator+ (const bbox& b, const bboxset& s) { - return bboxset(b) + s; -} - -// cost: O(1) -template -inline bboxset operator- (const bbox& b1, const bbox& b2) { - return bboxset::minus(b1,b2); -} - -// cost: O(n^2) -template -inline bboxset operator- (const bbox& b, const bboxset& s) { - return bboxset::minus(b,s); -} - - - -template -inline bboxset operator| (const bbox& b, const bboxset& s) { - return s | b; -} - -template -inline bboxset operator& (const bbox& b, const bboxset& s) { - return s & b; -} - - - -template -inline bool operator== (const bbox& b, const bboxset& s) -{ - return bboxset(b) == s; -} - -template -inline bool operator!= (const bbox& b, const bboxset& s) -{ - return bboxset(b) != s; -} - -template -inline bool operator< (const bbox& b, const bboxset& s) -{ - return bboxset(b) < s; -} - -template -inline bool operator<= (const bbox& b, const bboxset& s) -{ - return bboxset(b) <= s; -} - -template -inline bool operator> (const bbox& b, const bboxset& s) -{ - return bboxset(b) > s; -} - -template -inline bool operator>= (const bbox& b, const bboxset& s) -{ - return bboxset(b) >= s; -} - - - -template -inline bool operator== (const bboxset& s, const bbox& b) -{ - return s == bboxset(b); -} - -template -inline bool operator!= (const bboxset& s, const bbox& b) -{ - return s != bboxset(b); -} - -template -inline bool operator< (const bboxset& s, const bbox& b) -{ - return s < bboxset(b); -} - -template -inline bool operator<= (const bboxset& s, const bbox& b) -{ - return s <= bboxset(b); -} - -template -inline bool operator> (const bboxset& s, const bbox& b) -{ - return s > bboxset(b); -} - -template -inline bool operator>= (const bboxset& s, const bbox& b) -{ - return s >= bboxset(b); -} - - - -// Memory usage -template -inline size_t memoryof (bboxset const & s) -{ return s.memory(); } - - - -// Input -template -inline istream& operator>> (istream& is, bboxset& s) { - return s.input(is); -} - - - -// Output -template -inline ostream& operator<< (ostream& os, const bboxset& 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 +#include +#include +#include +#include +#include +#include + +#include "defs.hh" + +#include "bboxset1.hh" + +using namespace std; + + + +namespace bboxset1 { + + + +template class bboxset; +template void bboxset::serialise (set >& s) const; +template void bboxset::serialise (vector >& v) const; +template size_t memoryof (const bboxset& s); +template istream& operator>> (istream& is, bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& s); + +template class bboxset; +template void bboxset::serialise (set >& s) const; +template void bboxset::serialise (vector >& v) const; +template size_t memoryof (const bboxset& s); +template istream& operator>> (istream& is, bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& s); + +template class bboxset; +template void bboxset::serialise (set >& s) const; +template void bboxset::serialise (vector >& f) const; +template size_t memoryof (const bboxset& s); +template istream& operator>> (istream& is, bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& s); + +} // namespace bboxset1 + + + +#include "dh.hh" +#include "region.hh" + +namespace bboxset1 { + +template bboxset::bboxset (const vector& vb, const bbox dh::full_dboxes::* const v); +template bboxset::bboxset (const vector& vb, const bboxset dh::full_dboxes::* const v); +template bboxset::bboxset (const vector& vb, const bbox 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 +#include +#include +#include +#include + +#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 class bboxset; + +// template +// bboxset operator+ (const bbox& b1, const bbox& b2); +// template +// bboxset operator+ (const bbox& b, const bboxset& s); + +// template +// bboxset operator- (const bbox& b1, const bbox& b2); +// template +// bboxset operator- (const bbox& b, const bboxset& s); + +// Input +template +istream& operator>> (istream& is, bboxset& s); + +// Output +template +ostream& operator<< (ostream& os, const bboxset& s); + + + +// Bounding box set class +template +class bboxset { + + // Cost annotations depend on the number of bset elements n, and + // assume that normalization is skipped. + + struct skip_normalize_t { + bboxset& s; + bool const saved_skip_normalize; + skip_normalize_t(bboxset& 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 box; +#ifdef BBOXSET_SET + //S typedef set bset; +#endif +#ifdef BBOXSET_LIST + typedef list bset; +#endif +#ifdef BBOXSET_VECTOR + typedef vector 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& lb); + bboxset (const set& sb); + bboxset (const vector& vb); + bboxset (const vector >& vlb); + template + bboxset (const vector& vb, const bbox U::* const v); + template + bboxset (const vector& 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& 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- (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- (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& lo, const vect& hi) + const CCTK_MEMBER_ATTRIBUTE_PURE; + bboxset expand (const vect,2>& lohi) const + { return expand (lohi[0], lohi[1]); } + + /** Shift the bboxset by multiples of the stride. */ + bboxset shift (const vect& 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& lo, const vect& hi, + const vect& denom) const CCTK_MEMBER_ATTRIBUTE_PURE; + // cost: O(n^2) in general, but only O(n) for shifting + bboxset expand (const vect,2>& lohi, const vect& 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& v, const vect& 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 + 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 +inline bboxset operator+ (const bbox& b1, const bbox& b2) { + return bboxset(b1) + bboxset(b2); +} + +template +inline bboxset operator+ (const bbox& b, const bboxset& s) { + return bboxset(b) + s; +} + +// cost: O(1) +template +inline bboxset operator- (const bbox& b1, const bbox& b2) { + return bboxset::minus(b1,b2); +} + +// cost: O(n^2) +template +inline bboxset operator- (const bbox& b, const bboxset& s) { + return bboxset::minus(b,s); +} + + + +template +inline bboxset operator| (const bbox& b, const bboxset& s) { + return s | b; +} + +template +inline bboxset operator& (const bbox& b, const bboxset& s) { + return s & b; +} + + + +template +inline bool operator== (const bbox& b, const bboxset& s) +{ + return bboxset(b) == s; +} + +template +inline bool operator!= (const bbox& b, const bboxset& s) +{ + return bboxset(b) != s; +} + +template +inline bool operator< (const bbox& b, const bboxset& s) +{ + return bboxset(b) < s; +} + +template +inline bool operator<= (const bbox& b, const bboxset& s) +{ + return bboxset(b) <= s; +} + +template +inline bool operator> (const bbox& b, const bboxset& s) +{ + return bboxset(b) > s; +} + +template +inline bool operator>= (const bbox& b, const bboxset& s) +{ + return bboxset(b) >= s; +} + + + +template +inline bool operator== (const bboxset& s, const bbox& b) +{ + return s == bboxset(b); +} + +template +inline bool operator!= (const bboxset& s, const bbox& b) +{ + return s != bboxset(b); +} + +template +inline bool operator< (const bboxset& s, const bbox& b) +{ + return s < bboxset(b); +} + +template +inline bool operator<= (const bboxset& s, const bbox& b) +{ + return s <= bboxset(b); +} + +template +inline bool operator> (const bboxset& s, const bbox& b) +{ + return s > bboxset(b); +} + +template +inline bool operator>= (const bboxset& s, const bbox& b) +{ + return s >= bboxset(b); +} + + + +// Memory usage +template +inline size_t memoryof (bboxset const & s) +{ return s.memory(); } + + + +// Input +template +inline istream& operator>> (istream& is, bboxset& s) { + return s.input(is); +} + + + +// Output +template +inline ostream& operator<< (ostream& os, const bboxset& 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 +bboxset::bboxset () + : skip_normalize(false) +{ + assert (invariant()); +} + +template +bboxset::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 +bboxset::bboxset (const bboxset& s) + : bs(s.bs), skip_normalize(false) +{ + assert (invariant()); +} + +template +bboxset::bboxset (const list& lb) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename list::const_iterator + li = lb.begin(), le = lb.end(); li != le; ++ li) + { + *this |= *li; + } +} + +template +bboxset::bboxset (const set& sb) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename set::const_iterator + vi = sb.begin(), ve = sb.end(); vi != ve; ++ vi) + { + *this |= *vi; + } +} + +template +bboxset::bboxset (const vector& vb) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename vector::const_iterator + vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi) + { + *this |= *vi; + } +} + +template +bboxset::bboxset (const vector >& vlb) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename vector >::const_iterator + vli = vlb.begin(), vle = vlb.end(); vli != vle; ++ vli) + { + *this |= bboxset (*vli); + } +} + +template +template +bboxset::bboxset (const vector& vb, const bbox U::* const v) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename vector::const_iterator + vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi) + { + *this |= (*vi).*v; + } +} + +template +template +bboxset::bboxset (const vector& vb, const bboxset U::* const v) + : skip_normalize(false) +{ + SKIP_NORMALIZE(*this); + for (typename vector::const_iterator + vi = vb.begin(), ve = vb.end(); vi != ve; ++ vi) + { + *this |= (*vi).*v; + } +} + +template +bboxset bboxset::poison () { + return bboxset (bbox::poison()); +} + + + +// Invariant +template +bool bboxset::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 +void bboxset::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 buf; + typedef vector 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; dsize(); + + // Can't use operators on *this since these would call normalize again + // assert (*this == oldbs); + assert (newsize == oldsize); +} + + + +// Accessors +// cost: O(n) +template +typename bboxset::size_type bboxset::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::max() - bsz >= s); + s += bsz; + } + return s; +} + + + +// Queries + +// Containment +// cost: O(n) +template +bool bboxset::contains (const vect& 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 +bool bboxset::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 +bboxset& bboxset::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 +bboxset& bboxset::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 +bboxset& bboxset::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 +bboxset bboxset::operator+ (const box& b) const +{ + bboxset r(*this); + r += b; + assert (r.invariant()); + return r; +} + +// cost: O(n) +template +bboxset bboxset::operator+ (const bboxset& s) const +{ + bboxset r(*this); + r += s; + assert (r.invariant()); + return r; +} + + + +// Union +// cost: O(n) +template +bboxset& bboxset::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 +bboxset& bboxset::operator|= (const bboxset& s) +{ + bboxset tmp = s - *this; + add_transfer (tmp); + assert (invariant()); + return *this; +} + +// cost: O(n) +template +bboxset bboxset::operator| (const box& b) const +{ + bboxset r(*this); + r |= b; + assert (r.invariant()); + return r; +} + +// cost: O(n^2) +template +bboxset bboxset::operator| (const bboxset& s) const +{ + bboxset r(*this); + r |= s; + assert (r.invariant()); + return r; +} + + + +// Intersection +// cost: O(n) +template +bboxset bboxset::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 +bboxset bboxset::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 +bboxset& bboxset::operator&= (const box& b) +{ + *this = *this & b; + assert (invariant()); + return *this; +} + +// cost: O(n) +template +bboxset& bboxset::operator&= (const bboxset& s) +{ + *this = *this & s; + assert (invariant()); + return *this; +} + + + +// Difference +// cost: O(1) +template +bboxset bboxset::minus (const bbox& b1, const bbox& b2) +{ + assert (b1.is_aligned_with(b2)); + if (b1.empty()) return bboxset(); + if (not b1.intersects(b2)) return bboxset(b1); + vect const b2lo = max (b1.lower(), b2.lower()); + vect const b2up = min (b1.upper(), b2.upper()); + vect const & b1lo = b1.lower(); + vect const & b1up = b1.upper(); + vect const & str = b1.stride(); + bboxset r; + SKIP_NORMALIZE(r); + { + for (int d=0; d lb, ub; + bbox b; + for (int dd=0; ddd) { + lb[dd] = b1lo[dd]; + ub[dd] = b1up[dd]; + } + } + lb[d] = b1lo[d]; + ub[d] = b2lo[d] - str[d]; + b = bbox(lb,ub,str); + r += b; + lb[d] = b2up[d] + str[d]; + ub[d] = b1up[d]; + b = bbox(lb,ub,str); + r += b; + } + } + assert (r.invariant()); + return r; +} + +// cost: O(n) +template +bboxset bboxset::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 +bboxset& bboxset::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 +bboxset bboxset::operator- (const bboxset& s) const +{ + bboxset r(*this); + r -= s; + assert (r.invariant()); + return r; +} + +// cost: O(n^2) +template +bboxset bboxset::minus (const bbox& b, const bboxset& s) +{ + bboxset r = bboxset(b) - s; + assert (r.invariant()); + return r; +} + + + +// cost: O(n) +template +typename bboxset::box bboxset::container () const +{ + box b; + for (const_iterator bi=begin(), be=end(); bi!=be; ++bi) { + b = b.expanded_containing(*bi); + } + return b; +} + +template +bboxset bboxset::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 +bboxset bboxset::expand (const vect& lo, const vect& 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 +bboxset bboxset::expand (const vect& lo, const vect& hi, + const vect& denom) const +{ + assert (all(denom > vect(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 +bboxset bboxset::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 +bboxset bboxset::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 +bool bboxset::operator<= (const bboxset& s) const +{ + return (*this - s).empty(); +} + +template +bool bboxset::operator< (const bboxset& s) const +{ + return (*this - s).empty() && not (s - *this).empty(); +} + +template +bool bboxset::operator>= (const bboxset& s) const +{ + return s <= *this; +} + +template +bool bboxset::operator> (const bboxset& s) const +{ + return s < *this; +} + +template +bool bboxset::operator== (const bboxset& s) const +{ + return (*this <= s) && (*this >= s); +} + +template +bool bboxset::operator!= (const bboxset& s) const +{ + return not (*this == s); +} + + + +// Serialise +template +template +void bboxset::serialise (C& out) const +{ + for (bboxset::const_iterator it = begin(); it != end(); ++it) { + out.insert(out.end(), *it); + } +} + + + +// Input +template +istream& bboxset::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 +ostream& bboxset::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 + +#ifndef CARPET_NO_BBOXSET2 + +#include "bboxset2.hh" + + + +namespace bboxset2 { + + template class bboxset; + + template class bboxset; + template void bboxset::serialise(set&) const; + + template class bboxset; + template void bboxset::serialise(set&) const; + + template class bboxset; + template void bboxset::serialise(set&) const; + template void bboxset::serialise(vector&) const; + +} // namespace bboxset2 + + + +#include "dh.hh" +#include "region.hh" + +namespace bboxset2 { + + template bboxset::bboxset(const vector&, const bbox dh::full_dboxes::*); + template bboxset::bboxset(const vector&, const bboxset dh::full_dboxes::*); + template bboxset::bboxset(const vector&, 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "bbox.hh" +#include "defs.hh" +#include "vect.hh" + +using namespace std; + + + +#ifndef CARPET_NO_BBOXSET2 + +namespace bboxset2 { + + + +template +class bboxset { + template friend class bboxset; + typedef ::vect vect; + typedef ::bbox bbox; + typedef ::vect vect1; + typedef ::bbox bbox1; + typedef bboxset bboxset1; + +#if 0 && !defined __ICC + template + using ptr = unique_ptr; +#else + template + using ptr = shared_ptr; +#endif + + typedef map> subsets_t; + subsets_t subsets; + + vect stride, offset; + + bool is_poison_; + + template + 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 + 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::max(); + const T next_pos1 = iter1!=end1 ? iter1->first : numeric_limits::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 + 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(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 + bboxset(const C& elts); + + /** Create set from container of structs containing a bbox, bboxset, + or container thereof */ + template + 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& 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 + void serialise(C& out) const; + + /** Iterate over a serialised set */ +private: + typedef vector 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 +class bboxset { + template friend class bboxset; + typedef ::vect vect; + typedef ::bbox 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& lohi) const + { + return expand(lohi[0], lohi[1]); + } + + template + 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 +inline bboxset operator&(const bboxset& bs, const bbox& b) +{ + return bs & bboxset(b); +} +template +inline bboxset operator&(const bbox& b, const bboxset& bs) +{ + return bboxset(b) & bs; +} +// Note: bbox & bbox -> bbox + +template +inline bboxset operator^(const bboxset& bs, const bbox& b) +{ + return bs ^ bboxset(b); +} +template +inline bboxset operator^(const bbox& b, const bboxset& bs) +{ + return bboxset(b) ^ bs; +} +template +inline bboxset operator^(const bbox& b1, const bbox& b2) +{ + return bboxset(b1) ^ bboxset(b2); +} + +template +inline bboxset operator|(const bboxset& bs, const bbox& b) +{ + return bs | bboxset(b); +} +template +inline bboxset operator|(const bbox& b, const bboxset& bs) +{ + return bboxset(b) | bs; +} +template +inline bboxset operator|(const bbox& b1, const bbox& b2) +{ + return bboxset(b1) | bboxset(b2); +} + +template +inline bboxset operator+(const bboxset& bs, const bbox& b) +{ + return bs + bboxset(b); +} +template +inline bboxset operator+(const bbox& b, const bboxset& bs) +{ + return bboxset(b) + bs; +} +template +inline bboxset operator+(const bbox& b1, const bbox& b2) +{ + return bboxset(b1) + bboxset(b2); +} + +template +inline bboxset operator-(const bboxset& bs, const bbox& b) +{ + return bs - bboxset(b); +} +template +inline bboxset operator-(const bbox& b, const bboxset& bs) +{ + return bboxset(b) - bs; +} +template +inline bboxset operator-(const bbox& b1, const bbox& b2) +{ + return bboxset(b1) - bboxset(b2); +} + +template +inline bool operator==(const bbox& b, const bboxset& bs) +{ + return bboxset(b) == bs; +} +template +inline bool operator==(const bboxset& bs, const bbox& b) +{ + return bs == bboxset(b); +} + +template +inline bool operator!=(const bbox& b, const bboxset& bs) +{ + return bboxset(b) != bs; +} +template +inline bool operator!=(const bboxset& bs, const bbox& b) +{ + return bs != bboxset(b); +} + +template +inline bool operator<=(const bbox& b, const bboxset& bs) +{ + return bboxset(b) <= bs; +} +template +inline bool operator<=(const bboxset& bs, const bbox& b) +{ + return bs <= bboxset(b); +} + +template +inline bool operator<(const bbox& b, const bboxset& bs) +{ + return bboxset(b) < bs; +} +template +inline bool operator<(const bboxset& bs, const bbox& b) +{ + return bs < bboxset(b); +} + +template +inline bool operator>=(const bbox& b, const bboxset& bs) +{ + return bboxset(b) >= bs; +} +template +inline bool operator>=(const bboxset& bs, const bbox& b) +{ + return bs >= bboxset(b); +} + +template +inline bool operator>(const bbox& b, const bboxset& bs) +{ + return bboxset(b) > bs; +} +template +inline bool operator>(const bboxset& bs, const bbox& b) +{ + return bs > bboxset(b); +} + + + +/** Memory usage */ +template +inline size_t memoryof(bboxset const& bs) +{ + return bs.memory(); +} + +/** Input */ +template +inline istream& operator>>(istream& is, bboxset& bs) +{ + return bs.input(is); +} + +/** Output */ +template +inline ostream& operator<<(ostream& os, const bboxset& bs) +{ + return bs.output(os); +} + + + + ////////////////////////////////////////////////////////////////////////////// + // bboxset + ////////////////////////////////////////////////////////////////////////////// + + /** Copy constructor */ + template + bboxset::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(*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 + bboxset& bboxset::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(*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 + bboxset::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(b1); + const auto hi_subsetp = make_shared(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 + template + bboxset::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 + template + bboxset::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 + vect bboxset::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 + T bboxset::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 + T bboxset::size() const + { + assert(not is_poison()); + T total_size = 0; // accumulated total number of points + T old_pos = numeric_limits::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 + bbox bboxset::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 + bool bboxset::operator==(const bboxset& other) const + { + return (*this ^ other).empty(); + } + + /** Test for is-subset-of */ + template + bool bboxset::operator<=(const bboxset& other) const + { + return (*this | other) == other; + } + + /** Test for is-strict-subset-of */ + template + bool bboxset::operator<(const bboxset& other) const + { + return *this != other and *this <= other; + } + + + + /** Symmetric set difference */ + template + bboxset bboxset::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 + bboxset& bboxset::operator^=(const bboxset& other) + { + bboxset res = *this; + return *this = *this ^ other; + } + + /** Set intersection */ + template + bboxset bboxset::operator&(const bboxset& other) const + { + return binary_operator + ([](const bboxset1& set0, const bboxset1& set1) { return set0 & set1; }, + other); + } + + /** Set intersection */ + template + bboxset& bboxset::operator&=(const bboxset& other) + { + return *this = *this & other; + } + + /** Set Union */ + template + bboxset bboxset::operator|(const bboxset& other) const + { + return binary_operator + ([](const bboxset1& set0, const bboxset1& set1) { return set0 | set1; }, + other); + } + + /** Set union */ + template + bboxset& bboxset::operator|=(const bboxset& other) + { + return *this = *this | other; + } + + /** Symmetric set union */ + template + bboxset bboxset::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 + bboxset& bboxset::operator+=(const bboxset& other) + { + return *this = *this + other; + } + + /** Set difference */ + template + bboxset bboxset::operator-(const bboxset& other) const + { + return binary_operator + ([](const bboxset1& set0, const bboxset1& set1) { return set0 - set1; }, + other); + } + + /** Set difference */ + template + bboxset& bboxset::operator-=(const bboxset& other) + { + return *this = *this - other; + } + + + + /** Shift all points */ + template + bboxset bboxset::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(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 + bboxset bboxset:: 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 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 + bool bboxset::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 + bboxset bboxset::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 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 + bboxset bboxset::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 + bboxset bboxset::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 bs; + serialise(bs); + bboxset res; + for (const auto& b: bs) { + res |= b.anti_contracted_for(target); + } + return res; + } + + + + /** Serialise the set */ + template + template + void bboxset::serialise(C& out) const + { + assert(not is_poison()); + typedef map subboxes_t; + typedef set 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_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 + size_t bboxset::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 + istream& bboxset::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:"); + set 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 + ostream& bboxset::output(ostream& os) const + { + assert(not is_poison()); + T Tdummy; + set bs; + serialise(bs); + return os << "bboxset<" << typestring(Tdummy) << "," << D << ">(" + << "set:" << bs << "," + << "stride:" << stride << "," + << "offset:" << offset << ")"; + } + + + + ////////////////////////////////////////////////////////////////////////////// + // bboxset + ////////////////////////////////////////////////////////////////////////////// + + template + istream& bboxset::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 + ostream& bboxset::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 #include #include +#include #include #include #include +#include #include "bbox.hh" #include "defs.hh" @@ -283,6 +285,14 @@ ostream& output (ostream& os, const set& s) { return os; } +#ifndef CARPET_NO_BBOXSET2 +// Shared pointer output +template +ostream& output (ostream& os, const shared_ptr& s) { + return os << "(&" << *s.get() << ")"; +} +#endif + // Stack output template ostream& output (ostream& os, const stack& s) { @@ -327,13 +337,11 @@ ostream& output (ostream& os, const vector& v) { template int ipow (int x, int y); template CCTK_REAL ipow (CCTK_REAL x, int y); -template vect ipow (vect x, int y); +//template vect ipow (vect x, int y); template vect ipow (vect x, int y); -template bool equals (vector const& v, vector const& w); - -template size_t memoryof (list const & l); -template size_t memoryof (list const & l); +//template size_t memoryof (list const & l); +//template size_t memoryof (list const & l); template size_t memoryof (list const & l); template size_t memoryof (list const & l); template size_t memoryof (list const & l); @@ -342,49 +350,61 @@ template size_t memoryof (list const & l); template size_t memoryof (stack const & s); template size_t memoryof (vector const & v); template size_t memoryof (vector const & v); -template size_t memoryof (vector const & v); -template size_t memoryof (vector const & v); -template size_t memoryof (vector const & v); +//template size_t memoryof (vector const & v); +template size_t memoryof (vector > const & v); +template size_t memoryof (vector > const & v); +template size_t memoryof (vector > const & v); +template size_t memoryof (vector > const & v); +#ifndef CARPET_NO_BBOXSET2 +template size_t memoryof (vector> const & v); +#endif template size_t memoryof (vector const & v); template size_t memoryof (vector const & v); -template size_t memoryof (vector *> const & f); -template size_t memoryof (vector const & v); -template size_t memoryof (vector const & v); +template size_t memoryof (vector*> const & f); +//template size_t memoryof (vector const & v); +//template size_t memoryof (vector const & v); template size_t memoryof (vector const & v); template size_t memoryof (vector > const & v); template size_t memoryof (vector > const & v); template size_t memoryof (vector > const & v); template size_t memoryof (vector > const & v); -template size_t memoryof (vector > const & v); +//template size_t memoryof (vector > const & v); template size_t memoryof (vector > const & v); -template size_t memoryof (vector > const & v); -template size_t memoryof (vector > const & v); -template size_t memoryof (vector > const & v); +//template size_t memoryof (vector > const & v); +//template size_t memoryof (vector > const & v); +//template size_t memoryof (vector > const & v); template size_t memoryof (vector > > const & v); -template size_t memoryof (vector > > const & v); -template size_t memoryof (vector > > const & v); -template size_t memoryof (vector > > const & v); +//template size_t memoryof (vector > > const & v); +//template size_t memoryof (vector > > const & v); +//template size_t memoryof (vector > > const & v); template size_t memoryof (vector > > const & v); template size_t memoryof (vector > > const & v); template size_t memoryof (vector > > const & v); -template size_t memoryof (vector > > const & v); +//template size_t memoryof (vector > > const & v); template size_t memoryof (vector > > > const & v); -template istream& input (istream& os, list& l); -template istream& input (istream& os, set& s); +//template istream& input (istream& os, list& l); +template istream& input (istream& os, set >& s); +template istream& input (istream& os, set >& s); +template istream& input (istream& os, set >& s); template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); +//template istream& input (istream& os, vector& v); +template istream& input (istream& os, vector >& v); +template istream& input (istream& os, vector >& v); +template istream& input (istream& os, vector >& v); +//template istream& input (istream& os, vector& v); +template istream& input (istream& os, vector >& v); +#ifndef CARPET_NO_BBOXSET2 +template istream& input (istream& os, vector>& v); +#endif template istream& input (istream& os, vector& v); template istream& input (istream& os, vector& v); template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); -template istream& input (istream& os, vector& v); +//template istream& input (istream& os, vector& v); +//template istream& input (istream& os, vector& v); template istream& input (istream& os, vector& v); template istream& input (istream& os, vector >& v); -template istream& input (istream& os, vector >& v); +//template istream& input (istream& os, vector >& v); template istream& input (istream& os, vector >& v); template istream& input (istream& os, vector >& v); template istream& input (istream& os, vector >& v); @@ -393,31 +413,53 @@ template istream& input (istream& os, vector >& v); template istream& input (istream& os, vector > >& v); template istream& input (istream& os, vector > >& v); -template ostream& output (ostream& os, const list& l); -template ostream& output (ostream& os, const list& l); -template ostream& output (ostream& os, const map& m); -template ostream& output (ostream& os, const set& s); -template ostream& output (ostream& os, const set& s); -template ostream& output (ostream& os, const stack& s); +//template ostream& output (ostream& os, const list& l); +//template ostream& output (ostream& os, const list& l); +#ifndef CARPET_NO_BBOXSET2 +//template ostream& output (ostream& os, const map >>& m); +//template ostream& output (ostream& os, const map >>& m); +//template ostream& output (ostream& os, const map >>& m); +//template ostream& output (ostream& os, const map >>& m); +#endif +//template ostream& output (ostream& os, const map& m); +//template ostream& output (ostream& os, const pair,int>& p); +//template ostream& output (ostream& os, const set >& s); +template ostream& output (ostream& os, const set >& s); +template ostream& output (ostream& os, const set >& s); +template ostream& output (ostream& os, const set >& s); +//template ostream& output (ostream& os, const set >& s); +#ifndef CARPET_NO_BBOXSET2 +//template ostream& output (ostream& os, const shared_ptr >& s); +#endif +//template ostream& output (ostream& os, const stack& s); template ostream& output (ostream& os, const vector& v); template ostream& output (ostream& os, const vector& v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector& v); +template ostream& output (ostream& os, const vector >& v); +template ostream& output (ostream& os, const vector >& v); +template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector& v); +template ostream& output (ostream& os, const vector >& v); +#ifndef CARPET_NO_BBOXSET2 +template ostream& output (ostream& os, const vector>& v); +#endif template ostream& output (ostream& os, const vector& v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector& v); +//template ostream& output (ostream& os, const vector& v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector & v); -template ostream& output (ostream& os, const vector & v); -template ostream& output (ostream& os, const vector & v); -template ostream& output (ostream& os, const vector & v); -template ostream& output (ostream& os, const vector & v); +//template ostream& output (ostream& os, const vector & v); +//template ostream& output (ostream& os, const vector & v); +//template ostream& output (ostream& os, const vector & v); +//template ostream& output (ostream& os, const vector & v); +//template ostream& output (ostream& os, const vector & v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector& v); +//template ostream& output (ostream& os, const vector& v); template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector >& v); +//template ostream& output (ostream& os, const vector >& v); +//template ostream& output (ostream& os, const vector,int> >& v); +//template ostream& output (ostream& os, const vector,int> >& v); +//template ostream& output (ostream& os, const vector,int> >& v); +//template ostream& output (ostream& os, const vector,int> >& v); template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector >& v); @@ -425,16 +467,16 @@ template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector > & b); -template ostream& output (ostream& os, const vector > & b); +//template ostream& output (ostream& os, const vector > & b); template ostream& output (ostream& os, const vector > & b); -template ostream& output (ostream& os, const vector > & b); -template ostream& output (ostream& os, const vector > & b); +//template ostream& output (ostream& os, const vector > & b); +//template ostream& output (ostream& os, const vector > & b); template ostream& output (ostream& os, const vector >& v); template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > > & b); -template ostream& output (ostream& os, const vector > > & b); -template ostream& output (ostream& os, const vector > > & b); +//template ostream& output (ostream& os, const vector > >& v); +//template ostream& output (ostream& os, const vector > > & b); +//template ostream& output (ostream& os, const vector > > & b); +//template ostream& output (ostream& os, const vector > > & b); template ostream& output (ostream& os, const vector > > & b); template ostream& output (ostream& os, const vector > > & b); template ostream& output (ostream& os, const vector > >& 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 #include #include +#include #include #include #include @@ -62,13 +63,16 @@ char const * const eol = "\n"; #else # define AT(index) operator[](index) #endif - + // Some shortcuts for type names template class vect; template class bbox; -template class bboxset; +namespace bboxset1 { template class bboxset; } +#ifndef CARPET_NO_BBOXSET2 +namespace bboxset2 { template class bboxset; } +#endif template class fulltree; typedef vect bvect; @@ -78,15 +82,18 @@ typedef vect rvect; typedef bbox ibbox; typedef bbox jbbox; typedef bbox rbbox; -typedef bboxset ibset; - +namespace bboxset1 { typedef bboxset ibset; } +#ifndef CARPET_NO_BBOXSET2 +namespace bboxset2 { typedef bboxset ibset; } +#endif + // (Try to replace these by b2vect and i2vect) -typedef vect,dim> bbvect; -typedef vect,dim> iivect; +typedef vect,dim> bbvect; +typedef vect,dim> iivect; typedef vect,dim> jjvect; -typedef vect,2> b2vect; -typedef vect,2> i2vect; +typedef vect,2> b2vect; +typedef vect,2> i2vect; typedef vect,2> j2vect; typedef vect,2> r2vect; @@ -331,6 +338,9 @@ template ostream& output (ostream& os, const list& l); template ostream& output (ostream& os, const map& m); template ostream& output (ostream& os, const pair& p); template ostream& output (ostream& os, const set& s); +#ifndef CARPET_NO_BBOXSET2 +template ostream& output (ostream& os, const shared_ptr& s); +#endif template ostream& output (ostream& os, const stack& s); template ostream& output (ostream& os, const vector& v); @@ -344,11 +354,23 @@ inline ostream& operator<< (ostream& os, const map& m) { return output(os,m); } +template +inline ostream& operator<< (ostream& os, const pair& s) { + return output(os,s); +} + template inline ostream& operator<< (ostream& os, const set& s) { return output(os,s); } +#ifndef CARPET_NO_BBOXSET2 +template +inline ostream& operator<< (ostream& os, const shared_ptr& s) { + return output(os,s); +} +#endif + template inline ostream& operator<< (ostream& os, const stack& 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" -- cgit v1.2.3