diff options
Diffstat (limited to 'Carpet/CarpetLib/src')
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.cc | 79 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.hh | 21 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/defs.cc | 5 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 132 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/vect.cc | 35 |
5 files changed, 177 insertions, 95 deletions
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<T,D>::bboxset (const vector<list<box> >& vlb) { } } -template<class T, int D> +template<typename T, int D> +template<typename U> +bboxset<T,D>::bboxset (const vector<U>& vb, const bbox<T,D> U::* const v) { + for (typename vector<U>::const_iterator + vi = vb.begin(); vi != vb.end(); ++ vi) + { + *this |= (*vi).*v; + } +} + +template<typename T, int D> +template<typename U> +bboxset<T,D>::bboxset (const vector<U>& vb, const bboxset U::* const v) { + for (typename vector<U>::const_iterator + vi = vb.begin(); vi != vb.end(); ++ vi) + { + *this |= (*vi).*v; + } +} + +template<typename T, int D> bboxset<T,D> bboxset<T,D>::poison () { return bboxset (bbox<T,D>::poison()); } @@ -427,6 +447,53 @@ bboxset<T,D> bboxset<T,D>::minus (const bbox<T,D>& b, const bboxset<T,D>& s) { +template<typename T, int D> +typename bboxset<T,D>::box bboxset<T,D>::container () const { + box b; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + b = b.expanded_containing(*bi); + } + return b; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::pseudo_inverse (const int n) const { + assert (not empty()); + box const c = container().expand(n,n); + return c - *this; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& 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<typename T, int D> +bboxset<T,D> bboxset<T,D>::expanded_for (const box& b) const { + bboxset res; + for (const_iterator bi=begin(); bi!=end(); ++bi) { + res |= (*bi).expanded_for(b); + } + return res; +} + +template<typename T, int D> +bboxset<T,D> bboxset<T,D>::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<typename T, int D> bool bboxset<T,D>::operator<= (const bboxset<T,D>& s) const { @@ -512,3 +579,13 @@ ostream& bboxset<T,D>::output (ostream& os) const { template class bboxset<int,dim>; +template size_t memoryof (const bboxset<int,dim>& s); +template istream& operator>> (istream& is, bboxset<int,dim>& s); +template ostream& operator<< (ostream& os, const bboxset<int,dim>& s); + +#include "dh.hh" +#include "region.hh" + +template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bbox<int,dim> dh::full_dboxes::* const v); +template bboxset<int,dim>::bboxset (const vector<dh::full_dboxes>& vb, const bboxset<int,dim> dh::full_dboxes::* const v); +template bboxset<int,dim>::bboxset (const vector<region_t>& vb, const bbox<int,dim> region_t::* const v); diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index 85eb2b2c8..5624ef4b1 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.hh @@ -63,6 +63,10 @@ public: bboxset (const list<box>& lb); bboxset (const vector<list<box> >& vlb); + template<typename U> + bboxset (const vector<U>& vb, const bbox<T,D> U::* const v); + template<typename U> + bboxset (const vector<U>& vb, const bboxset U::* const v); static bboxset poison (); @@ -133,6 +137,23 @@ public: // friend bboxset operator- <T,D>(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<T,D>& lo, const vect<T,D>& hi) const; + bboxset expand (const vect<vect<T,D>,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<int,dim> ipow (vect<int,dim> x, int y); template vect<CCTK_REAL,dim> ipow (vect<CCTK_REAL,dim> x, int y); +template bool equals (vector<ibset> const& v, vector<ibset> const& w); + template size_t memoryof (list<ibbox> const & l); template size_t memoryof (list<ivect> const & l); template size_t memoryof (list<dh*> const & l); @@ -342,6 +344,7 @@ template size_t memoryof (vector<bool> const & v); template size_t memoryof (vector<int> const & v); template size_t memoryof (vector<CCTK_REAL> const & v); template size_t memoryof (vector<ibbox> const & v); +template size_t memoryof (vector<ibset> const & v); template size_t memoryof (vector<ivect> const & v); template size_t memoryof (vector<i2vect> const & v); template size_t memoryof (vector<fulltree <int,dim,pseudoregion_t> *> const & f); @@ -367,6 +370,7 @@ template istream& input (istream& os, vector<int>& v); template istream& input (istream& os, vector<CCTK_REAL>& v); template istream& input (istream& os, vector<ibbox>& v); template istream& input (istream& os, vector<rbbox>& v); +template istream& input (istream& os, vector<ibset>& v); template istream& input (istream& os, vector<ivect>& v); template istream& input (istream& os, vector<bbvect>& v); template istream& input (istream& os, vector<i2vect>& v); @@ -394,6 +398,7 @@ template ostream& output (ostream& os, const vector<int>& v); template ostream& output (ostream& os, const vector<CCTK_REAL>& v); template ostream& output (ostream& os, const vector<ibbox>& v); template ostream& output (ostream& os, const vector<rbbox>& v); +template ostream& output (ostream& os, const vector<ibset>& v); template ostream& output (ostream& os, const vector<ivect>& v); template ostream& output (ostream& os, const vector<rvect>& v); template ostream& output (ostream& os, const vector<bbvect>& v); 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<vect<bool,2>,dim>::output (ostream& os) const; template void vect<vect<int,2>,dim>::output (ostream& os) const; template void vect<vect<bool,dim>,2>::output (ostream& os) const; template void vect<vect<int,dim>,2>::output (ostream& os) const; +template void vect<vect<CCTK_REAL,dim>,2>::output (ostream& os) const; + + + +// Instantiate for bboxset class + +#define DEFINE_FAKE_VECT_OPERATIONS(T,D) \ +template<> vect<T,D> vect<T,D>::dir (const int d) { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq () { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq (const int n) { assert(0); } \ +template<> vect<T,D> vect<T,D>::seq (const int n, const int s) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator*= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator*= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator/= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator/= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator%= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator%= (const T&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator^= (const vect<T,D>&) { assert(0); } \ +template<> vect<T,D>& vect<T,D>::operator^= (const T&) { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator+ () const { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator- () const { assert(0); } \ +template<> vect<T,D> vect<T,D>::operator~ () const { assert(0); } \ +template class vect<T,D>; \ +template size_t memoryof (const vect<T,D>&); \ +template istream& operator>> (istream& is, vect<T,D>&); \ +template ostream& operator<< (ostream& os, const vect<T,D>&); + +typedef bboxset<int,dim> T1; +typedef vect<bboxset<int,dim>,2> T2; + +DEFINE_FAKE_VECT_OPERATIONS(T1,dim) +DEFINE_FAKE_VECT_OPERATIONS(T2,dim) + +#undef DEFINE_FAKE_VECT_OPERATIONS |