diff options
Diffstat (limited to 'Carpet/CarpetLib')
-rw-r--r-- | Carpet/CarpetLib/src/bbox.cc | 52 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bbox.hh | 46 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.cc | 114 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.hh | 29 |
4 files changed, 141 insertions, 100 deletions
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc index e26daf9c7..488b33598 100644 --- a/Carpet/CarpetLib/src/bbox.cc +++ b/Carpet/CarpetLib/src/bbox.cc @@ -13,42 +13,13 @@ using namespace std; -// Constructors -template<class T, int D> -bbox<T,D>::bbox (): _lower(T(1)), _upper(T(0)), _stride(T(1)) { } - -template<class T, int D> -bbox<T,D>::bbox (const bbox& b) - : _lower(b._lower), _upper(b._upper), _stride(b._stride) -{ } - -template<class T, int D> -bbox<T,D>& bbox<T,D>::operator= (const bbox& b) { - _lower=b._lower; _upper=b._upper; _stride=b._stride; - return *this; -} - -template<class T, int D> -bbox<T,D>::bbox (const vect<T,D>& lower_, const vect<T,D>& upper_, - const vect<T,D>& stride_) - : _lower(lower_), _upper(upper_), _stride(stride_) -{ - assert (all(_stride>T(0))); - assert (all((_upper-_lower)%_stride == T(0))); - if (numeric_limits<T>::is_integer && numeric_limits<T>::is_signed) { - // prevent accidental wrap-around - assert (all(_lower < numeric_limits<T>::max() / 2)); - assert (all(_lower > numeric_limits<T>::min() / 2)); - assert (all(_upper < numeric_limits<T>::max() / 2)); - assert (all(_upper > numeric_limits<T>::min() / 2)); - } -} - // Accessors template<class T, int D> T bbox<T,D>::size () const { if (empty()) return 0; -// return prod((shape()+stride()-1)/stride()); +#ifdef NDEBUG + return prod(shape()/stride()+1); +#else const vect<T,D> sh((shape()+stride()-T(1))/stride()); T sz = 1, max = numeric_limits<T>::max(); for (int d=0; d<D; ++d) { @@ -57,6 +28,7 @@ T bbox<T,D>::size () const { max /= sh[d]; } return sz; +#endif } // Queries @@ -114,15 +86,6 @@ bool bbox<T,D>::operator>= (const bbox& b) const { return b <= *this; } -// Intersection -template<class T, int D> -bbox<T,D> bbox<T,D>::operator& (const bbox& b) const { - assert (all(stride()==b.stride())); - vect<T,D> lo = max(lower(),b.lower()); - vect<T,D> up = min(upper(),b.upper()); - return bbox(lo,up,stride()); -} - // Containment template<class T, int D> bool bbox<T,D>::is_contained_in (const bbox& b) const { @@ -131,6 +94,13 @@ bool bbox<T,D>::is_contained_in (const bbox& b) const { return all(lower()>=b.lower() && upper()<=b.upper()); } +// Intersection +template<class T, int D> +bool bbox<T,D>::intersects (const bbox& b) const { + // no alignment check + return all(upper()>=b.lower() && lower()<=b.upper()); +} + // Alignment check template<class T, int D> bool bbox<T,D>::is_aligned_with (const bbox& b) const { diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index 8d0f8f54e..92563651b 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.hh @@ -1,7 +1,9 @@ #ifndef BBOX_HH #define BBOX_HH +#include <cassert> #include <iostream> +#include <limits> #include "defs.hh" #include "vect.hh" @@ -38,17 +40,40 @@ public: // Constructors /** Construct an empty bbox. */ - bbox (); + bbox () + : _lower(T(1)), _upper(T(0)), _stride(T(1)) + { + } /** Copy constructor. */ - bbox (const bbox& b); + bbox (const bbox& b) + : _lower(b._lower), _upper(b._upper), _stride(b._stride) + { + } /** Assignment operator. */ - bbox& operator= (const bbox& b); + bbox& operator= (const bbox& b) + { + _lower=b._lower; _upper=b._upper; _stride=b._stride; + return *this; + } /** Create a bbox from bounds and stride. */ - bbox (const vect<T,D>& lower, const vect<T,D>& upper, - const vect<T,D>& stride); + bbox (const vect<T,D>& lower_, + const vect<T,D>& upper_, + const vect<T,D>& stride_) + : _lower(lower_), _upper(upper_), _stride(stride_) + { + assert (all(_stride>T(0))); + assert (all((_upper-_lower)%_stride == T(0))); + if (numeric_limits<T>::is_integer and numeric_limits<T>::is_signed) { + // prevent accidental wrap-around + assert (all(_lower < numeric_limits<T>::max() / 2)); + assert (all(_lower > numeric_limits<T>::min() / 2)); + assert (all(_upper < numeric_limits<T>::max() / 2)); + assert (all(_upper > numeric_limits<T>::min() / 2)); + } +} // Accessors // (Don't return references; *this might be a temporary) @@ -88,11 +113,20 @@ public: /** Calculate the intersection (the set of common points) with the bbox b. */ - bbox operator& (const bbox& b) const; + bbox operator& (const bbox& b) const + { + assert (all(stride()==b.stride())); + vect<T,D> lo = max(lower(),b.lower()); + vect<T,D> up = min(upper(),b.upper()); + return bbox(lo,up,stride()); + } /** Find out whether this bbox is contained in the bbox b. */ bool is_contained_in (const bbox& b) const; + /** Find out whether this bbox intersects the bbox b. */ + bool intersects (const bbox& b) const; + /** Find out whether this bbox is aligned with the bbox b. ("aligned" means that both bboxes have the same stride and that their boundaries are commesurate.) */ diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 7725e899b..c36e8a524 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -21,7 +21,8 @@ bboxset<T,D>::bboxset () { template<class T, int D> bboxset<T,D>::bboxset (const box& b) { - if (!b.empty()) bs.insert(b); + //S if (!b.empty()) bs.insert(b); + if (!b.empty()) bs.push_back(b); assert (invariant()); } @@ -41,13 +42,13 @@ bboxset<T,D>::bboxset (const bset& bs_): bs(bs_) { template<class T, int D> bool bboxset<T,D>::invariant () const { // This is very slow when there are many bboxes -#if 0 +#if 0 && ! defined(CARPET_OPTIMISE) for (const_iterator bi=begin(); bi!=end(); ++bi) { if ((*bi).empty()) return false; - if (! (*bi).is_aligned_with(*bs.begin())) 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) & (*bi)).empty()) return false; + if ((*bi2).intersects(*bi)) return false; } } #endif @@ -69,15 +70,23 @@ void bboxset<T,D>::normalize () // boundaries aligned. for (int d=0; d<D; ++d) { // Find all boundaries - set<T> sbnds; + //S typedef set<T> buf; + typedef vector<T> buf; + buf sbnds; + sbnds.reserve (2 * bs.size()); for (typename bset::const_iterator si = bs.begin(); si != bs.end(); ++ si) { box const & b = * si; int const bstr = b.stride()[d]; int const blo = b.lower()[d]; int const bhi = b.upper()[d] + bstr; - sbnds.insert (blo); - sbnds.insert (bhi); + //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(); si != bs.end(); ++ si) { @@ -85,24 +94,30 @@ void bboxset<T,D>::normalize () int const bstr = b.stride()[d]; int const blo = b.lower()[d]; int const bhi = b.upper()[d] + bstr; - typename set<T>::const_iterator const ilo - = find (sbnds.begin(), sbnds.end(), blo); - typename set<T>::const_iterator const ihi - = find (sbnds.begin(), sbnds.end(), bhi); +// typename buf::const_iterator const ilo +// = find (sbnds.begin(), sbnds.end(), blo); +// typename buf::const_iterator const ihi +// = find (sbnds.begin(), sbnds.end(), bhi); + 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 set<T>::const_iterator curr = ilo; curr != ihi; ++ curr) { - typename set<T>::const_iterator next = curr; + 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 nb (b.lower().replace(d, nblo), b.upper().replace(d, nbhi - bstr), b.stride()); - nbs.insert (nb); + //S nbs.insert (nb); + nbs.push_back (nb); } } // Replace old set @@ -142,7 +157,8 @@ void bboxset<T,D>::normalize () b.stride()); bs.erase (si); nbs.erase (nsi); - bs.insert (cb); + //S bs.insert (cb); + bs.push_back (cb); goto done; } else if (bhi == nblo) { // Combine boxes, b < nb @@ -151,13 +167,15 @@ void bboxset<T,D>::normalize () b.stride()); bs.erase (si); nbs.erase (nsi); - bs.insert (cb); + //S bs.insert (cb); + bs.push_back (cb); goto done; } } } bs.erase (si); - nbs.insert (b); + //S nbs.insert (b); + nbs.push_back (b); done:; } bs.swap (nbs); @@ -188,22 +206,17 @@ T bboxset<T,D>::size () const { // Add (bboxes that don't overlap) template<class T, int D> -bboxset<T,D>& bboxset<T,D>::operator+= (const box& b) { - if (b.empty()) return *this; - // check for overlap - for (const_iterator bi=begin(); bi!=end(); ++bi) { - assert ((*bi & b).empty()); +bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) { + for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) { + *this += *bi; } - bs.insert(b); assert (invariant()); return *this; } template<class T, int D> -bboxset<T,D>& bboxset<T,D>::operator+= (const bboxset& s) { - for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) { - *this += *bi; - } +bboxset<T,D>& bboxset<T,D>::add_transfer (bboxset& s) { + bs.splice (bs.end(), s.bs); assert (invariant()); return *this; } @@ -239,14 +252,16 @@ bboxset<T,D> bboxset<T,D>::plus (const bbox<T,D>& b, const bboxset<T,D>& s) { // Union template<class T, int D> bboxset<T,D>& bboxset<T,D>::operator|= (const box& b) { - *this += b - *this; + bboxset tmp = b - *this; + add_transfer (tmp); assert (invariant()); return *this; } template<class T, int D> bboxset<T,D>& bboxset<T,D>::operator|= (const bboxset& s) { - *this += s - *this; + bboxset tmp = s - *this; + add_transfer (tmp); assert (invariant()); return *this; } @@ -290,7 +305,8 @@ bboxset<T,D> bboxset<T,D>::operator& (const bboxset& s) const { // walk all the bboxes for (const_iterator bi=s.begin(); bi!=s.end(); ++bi) { // insert the intersection with this bbox - r += *this & *bi; + bboxset tmp = *this & *bi; + r.add_transfer (tmp); } assert (r.invariant()); return r; @@ -317,8 +333,12 @@ template<class T, int D> bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) { assert (b1.is_aligned_with(b2)); if (b1.empty()) return bboxset<T,D>(); - if (b2.empty()) return bboxset<T,D>(b1); - const vect<T,D> str = b1.stride(); + if (not b1.intersects(b2)) return bboxset<T,D>(b1); + vect<T,D> const b2lo = max (b1.lower(), b2.lower()); + vect<T,D> const b2up = min (b1.upper(), b2.upper()); + vect<T,D> const & b1lo = b1.lower(); + vect<T,D> const & b1up = b1.upper(); + vect<T,D> const & str = b1.stride(); bboxset<T,D> r; for (int d=0; d<D; ++d) { // make resulting bboxes as large as possible in x-direction (for @@ -327,20 +347,20 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b1, const bbox<T,D>& b2) { bbox<T,D> b; for (int dd=0; dd<D; ++dd) { if (dd<d) { - lb[dd] = b2.lower()[dd]; - ub[dd] = b2.upper()[dd]; + lb[dd] = b2lo[dd]; + ub[dd] = b2up[dd]; } else if (dd>d) { - lb[dd] = b1.lower()[dd]; - ub[dd] = b1.upper()[dd]; + lb[dd] = b1lo[dd]; + ub[dd] = b1up[dd]; } } - lb[d] = b1.lower()[d]; - ub[d] = b2.lower()[d] - str[d]; - b = bbox<T,D>(lb,ub,str) & b1; + lb[d] = b1lo[d]; + ub[d] = b2lo[d] - str[d]; + b = bbox<T,D>(lb,ub,str); r += b; - lb[d] = b2.upper()[d] + str[d]; - ub[d] = b1.upper()[d]; - b = bbox<T,D>(lb,ub,str) & b1; + lb[d] = b2up[d] + str[d]; + ub[d] = b1up[d]; + b = bbox<T,D>(lb,ub,str); r += b; } assert (r.invariant()); @@ -354,18 +374,12 @@ bboxset<T,D> bboxset<T,D>::operator- (const box& b) const { // walk all my elements for (const_iterator bi=begin(); bi!=end(); ++bi) { // insert the difference with the bbox - r += *bi - b; + bboxset tmp = *bi - b; + r.add_transfer (tmp); } assert (r.invariant()); return r; } - -template<class T, int D> -bboxset<T,D>& bboxset<T,D>::operator-= (const box& b) { - *this = *this - b; - assert (invariant()); - return *this; -} template<class T, int D> bboxset<T,D>& bboxset<T,D>::operator-= (const bboxset& s) { diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index 44eee3c0a..5ad2dea94 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.hh @@ -3,6 +3,7 @@ #include <cassert> #include <iostream> +#include <list> #include <set> #include "bbox.hh" @@ -38,7 +39,8 @@ class bboxset { // Types typedef bbox<T,D> box; - typedef set<box> bset; + //S typedef set<box> bset; + typedef list<box> bset; // Fields bset bs; @@ -67,8 +69,24 @@ public: int setsize () const { return bs.size(); } // Add (bboxes that don't overlap) - bboxset& operator+= (const box& b); + bboxset& operator+= (const box& b) + { + if (b.empty()) return *this; + // This is very slow when there are many bboxes +#if 0 && ! defined(CARPET_OPTIMISE) + // check for overlap + for (const_iterator bi=begin(); bi!=end(); ++bi) { + assert (not (*bi).intersects(b)); + } +#endif + //S bs.insert(b); + bs.push_back(b); + assert (invariant()); + return *this; + } + bboxset& operator+= (const bboxset& s); + bboxset& add_transfer (bboxset& s); bboxset operator+ (const box& b) const; bboxset operator+ (const bboxset& s) const; static bboxset plus (const box& b1, const box& b2); @@ -90,7 +108,12 @@ public: // friend bboxset operator- <T,D>(const box& b1, const box& b2); static bboxset minus (const box& b1, const box& b2); bboxset operator- (const box& b) const; - bboxset& operator-= (const box& b); + bboxset& operator-= (const box& b) + { + *this = *this - b; + assert (invariant()); + return *this; + } bboxset& operator-= (const bboxset& s); bboxset operator- (const bboxset& s) const; // friend bboxset operator- <T,D>(const box& b, const bboxset& s); |