aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetLib/src
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetLib/src')
-rw-r--r--Carpet/CarpetLib/src/bbox.cc52
-rw-r--r--Carpet/CarpetLib/src/bbox.hh46
-rw-r--r--Carpet/CarpetLib/src/bboxset.cc114
-rw-r--r--Carpet/CarpetLib/src/bboxset.hh29
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-01 05:00:16 +0000