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