SpatialOps
MemoryWindow.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2014 The University of Utah
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to
6  * deal in the Software without restriction, including without limitation the
7  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8  * sell copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20  * IN THE SOFTWARE.
21  */
22 
23 #ifndef SpatialOps_MemoryWindow_h
24 #define SpatialOps_MemoryWindow_h
25 
26 #include <cstddef>
27 #include <vector>
28 #include <iterator>
29 # include <sstream>
30 # include <iostream>
31 
32 #include <spatialops/SpatialOpsConfigure.h>
33 
37 
38 #ifndef NDEBUG
39 # include <cassert>
40 # include <stdexcept>
41 #endif
42 
47 namespace SpatialOps{
48 
63  class MemoryWindow{
64 
65  friend std::ostream& operator<<( std::ostream&, const MemoryWindow& );
66 
67  IntVec nptsGlob_;
68  IntVec offset_;
69  IntVec extent_;
70 
71  public:
72 
79  MemoryWindow( const int npts[3],
80  const int offset[3],
81  const int extent[3] );
82 
89  MemoryWindow( const IntVec& npts,
90  const IntVec& offset,
91  const IntVec& extent );
92 
97  MemoryWindow( const int npts[3] );
98 
103  MemoryWindow( const IntVec& npts );
104 
105  // inline to improve performance since this is frequently called from nebo
106  MemoryWindow( const MemoryWindow& other )
107  : nptsGlob_( other.nptsGlob_ ),
108  offset_ ( other.offset_ ),
109  extent_ ( other.extent_ )
110  {}
111 
112  MemoryWindow& operator=( const MemoryWindow& other );
113 
114  ~MemoryWindow();
115 
121  inline int flat_index( IntVec loc ) const{
122 # ifndef NDEBUG
123  if( extent_[0]>1 ) assert( loc[0] < nptsGlob_[0] );
124  if( extent_[1]>1 ) assert( loc[1] < nptsGlob_[1] );
125  if( extent_[2]>1 ) assert( loc[2] < nptsGlob_[2] );
126 # endif
127  loc[0] = nptsGlob_[0] > 1 ? loc[0]+offset_[0] : 0;
128  loc[1] = nptsGlob_[1] > 1 ? loc[1]+offset_[1] : 0;
129  loc[2] = nptsGlob_[2] > 1 ? loc[2]+offset_[2] : 0;
130  return ijk_to_flat(nptsGlob_,loc);
131  }
132 
138  inline IntVec ijk_index_from_global( const int loc ) const{
139  return flat_to_ijk(nptsGlob_,loc);
140  }
141 
147  inline IntVec ijk_index_from_local( const int loc ) const{
148  return flat_to_ijk(extent_,loc);
149  }
150 
155  inline size_t glob_npts() const{ return nptsGlob_[0] * nptsGlob_[1] * nptsGlob_[2]; }
156 
161  inline size_t local_npts() const{ return extent_[0] * extent_[1] * extent_[2]; }
162 
163  inline size_t glob_dim( const size_t i ) const{ assert(i<3); return size_t(nptsGlob_[i]); }
164  inline size_t offset ( const size_t i ) const{ assert(i<3); return size_t(offset_[i] ); }
165  inline size_t extent ( const size_t i ) const{ assert(i<3); return size_t(extent_[i] ); }
166 
167  inline int& offset( const size_t i ){ assert(i<3); return offset_[i]; }
168  inline int& extent( const size_t i ){ assert(i<3); return extent_[i]; }
169 
170  inline IntVec& extent(){ return extent_; }
171  inline IntVec& offset(){ return offset_; }
172 
173  inline const IntVec& extent () const{ return extent_; }
174  inline const IntVec& offset () const{ return offset_; }
175  inline const IntVec& glob_dim() const{ return nptsGlob_; }
176 
180  inline bool operator==( const MemoryWindow& w ) const{
181  return (nptsGlob_ == w.nptsGlob_) &&
182  (extent_ == w.extent_ ) &&
183  (offset_ == w.offset_ );
184  }
185 
186  inline bool operator!=( const MemoryWindow& w ) const{
187  return (nptsGlob_ != w.nptsGlob_) ||
188  (extent_ != w.extent_ ) ||
189  (offset_ != w.offset_ );
190  }
191 
198  MemoryWindow reset_ghosts( GhostData const & oldGhosts,
199  GhostData const & newGhosts ) const
200  {
201  IntVec oldTotal = oldGhosts.get_minus() + oldGhosts.get_plus();
202  IntVec newTotal = newGhosts.get_minus() + newGhosts.get_plus();
203 
204  return MemoryWindow( glob_dim(),
205  offset() + oldGhosts.get_minus() - newGhosts.get_minus(),
206  extent() - oldTotal + newTotal );
207  }
208 
214  MemoryWindow remove_ghosts( GhostData const & ghosts ) const {
215  return (*this).reset_ghosts( ghosts, GhostData(0,0,0,0,0,0) );
216  }
217 
227  bool fits_in( MemoryWindow const & outer ) const
228  {
229  return ( glob_dim() == outer.glob_dim() &&
230  offset() >= outer.offset() &&
231  extent() <= outer.extent() );
232  }
233 
247  bool fits_between( MemoryWindow const & inner,
248  MemoryWindow const & outer ) const
249  {
250  assert( inner.fits_in(outer) );
251  return ( inner.fits_in(*this) && (*this).fits_in(outer) );
252  }
253 
258  inline std::string print() const
259  {
260  std::stringstream s;
261  s << "Offset: " << offset_ << std::endl
262  << "Extent: " << extent_ << std::endl;
263  return s.str();
264  }
265 
269  bool sanity_check() const;
270 
271  };
272 
273  //============================================================================
274 
290  inline int get_dim_with_ghost( const int nNoGhost,
291  const int minusGhost,
292  const int plusGhost,
293  const int bc )
294  {
295  return ( nNoGhost > 1 ? ( nNoGhost + minusGhost + plusGhost + bc ) : 1 );
296  }
297 
298  //------------------------------------------------------------------
299 
314  inline get_window_with_ghost( const IntVec& localDim,
315  const GhostData& ghost,
316  const BoundaryCellInfo& bc )
317  {
318  return MemoryWindow( IntVec( get_dim_with_ghost( localDim[0], ghost.get_minus(0), ghost.get_plus(0), bc.has_extra(0) ),
319  get_dim_with_ghost( localDim[1], ghost.get_minus(1), ghost.get_plus(1), bc.has_extra(1) ),
320  get_dim_with_ghost( localDim[2], ghost.get_minus(2), ghost.get_plus(2), bc.has_extra(2) ) ) );
321  }
322 
323  //============================================================================
324 
325  template<typename FieldType>
326  class ConstFieldIterator; // forward
327 
333  template<typename FieldType>
334  class FieldIterator : public std::iterator<std::random_access_iterator_tag, typename FieldType::value_type> {
335  friend class ConstFieldIterator<FieldType>;
337  typedef typename FieldType::value_type AtomicType;
338 
339  public:
340  FieldIterator()
341  : current_(NULL),
342  count_(0),
343  xIndex_(0), yIndex_(0), zIndex_(0),
344  yStep_(0), zStep_(0),
345  xExtent_(0), yExtent_(0), zExtent_(0),
346  xyExtent_(0)
347  {}
348 
349  FieldIterator( AtomicType * field_values,
350  const MemoryWindow & w )
351  : current_(field_values +
352  w.offset(0) +
353  w.offset(1) * w.glob_dim(0) +
354  w.offset(2) * w.glob_dim(0) * w.glob_dim(1)),
355  count_(0),
356  xIndex_(0), yIndex_(0), zIndex_(0),
357  yStep_(w.glob_dim(0) - w.extent(0)),
358  zStep_((w.glob_dim(1) - w.extent(1)) * w.glob_dim(0)),
359  xExtent_(w.extent(0)),
360  yExtent_(w.extent(1)),
361  zExtent_(w.extent(2)),
362  xyExtent_(w.extent(0) * w.extent(1))
363  {}
364 
365  //mutable dereference
366  inline AtomicType & operator*() {
367 # ifndef NDEBUG
368  if( count_ != (xIndex_ +
369  yIndex_ * xExtent_ +
370  zIndex_ * xyExtent_) ){
371  std::ostringstream msg;
372  msg << __FILE__ << " : " << __LINE__ << std::endl
373  << "iterator's internal count is off";
374  throw std::runtime_error(msg.str());
375  }
376  if( xIndex_ >= xExtent_ ||
377  yIndex_ >= yExtent_ ||
378  zIndex_ >= zExtent_ ||
379  xIndex_ < 0 ||
380  yIndex_ < 0 ||
381  zIndex_ < 0 ){
382  std::ostringstream msg;
383  msg << __FILE__ << " : " << __LINE__ << std::endl
384  << "iterator is in an invalid state for dereference";
385  throw std::runtime_error(msg.str());
386  }
387 # endif
388  return *current_;
389  }
390 
391  //immutable dereference
392  inline AtomicType const & operator*() const {
393 # ifndef NDEBUG
394  if(count_ != (xIndex_ +
395  yIndex_ * xExtent_ +
396  zIndex_ * xyExtent_) ){
397  std::ostringstream msg;
398  msg << __FILE__ << " : " << __LINE__ << std::endl
399  << "iterator's internal count is off";
400  throw std::runtime_error(msg.str());
401  }
402  if( xIndex_ >= xExtent_ ||
403  yIndex_ >= yExtent_ ||
404  zIndex_ >= zExtent_ ||
405  xIndex_ < 0 ||
406  yIndex_ < 0 ||
407  zIndex_ < 0 ){
408  std::ostringstream msg;
409  msg << __FILE__ << " : " << __LINE__ << std::endl
410  << "iterator is in an invalid state for dereference";
411  throw std::runtime_error(msg.str());
412  }
413 # endif
414  return *current_;
415  }
416 
417  //increment
418  inline Self & operator++() {
419  current_++; //xStep
420  count_++;
421  xIndex_++;
422  if(xIndex_ == xExtent_) {
423  current_ += yStep_; //yStep
424  xIndex_ = 0;
425  yIndex_++;
426  if(yIndex_ == yExtent_) {
427  current_ += zStep_; //zStep
428  yIndex_ = 0;
429  zIndex_++;
430  }
431  }
432  return *this;
433  }
434  inline Self operator++(int) { Self result = *this; ++(*this); return result; }
435 
436  //decrement
437  inline Self & operator--() {
438  current_--; //xStep
439  count_--;
440  xIndex_--;
441  if( xIndex_ == -1 ){
442  current_ -= yStep_; //yStep
443  xIndex_ = xExtent_ - 1;
444  yIndex_--;
445  if( yIndex_ == -1 ){
446  current_ -= zStep_; //zStep
447  yIndex_ = yExtent_ - 1;
448  zIndex_--;
449  }
450  }
451  return *this;
452  }
453  inline Self operator--(int) { Self result = *this; --(*this); return result; }
454 
455  //compound assignment
456  inline Self & operator+=(int change) {
457  //small change (only changes xIndex_)
458  if( (change == 0) || //no change
459  (change > 0 && //positive change
460  change < xExtent_ - xIndex_) ||
461  (change < 0 && //negative change
462  -change < xIndex_) ){
463  current_ += change;
464  xIndex_ += change;
465  count_ += change;
466  }
467  //bigger change (changes yIndex_ and/or zIndex_)
468  else {
469  current_ += (change + //xStep
470  yStep_ * (((count_ + change) / xExtent_ ) - (count_ / xExtent_)) +
471  zStep_ * (((count_ + change) / xyExtent_) - (count_ /xyExtent_)));
472  count_ += change;
473  xIndex_ = count_ % xExtent_;
474  yIndex_ = (count_ % xyExtent_) / xExtent_;
475  zIndex_ = count_ / xyExtent_;
476  }
477  return *this;
478  }
479  inline Self & operator-=(int change) { return *this += -change; }
480 
481  //addition/subtraction
482  inline Self operator+ (int change) const { Self result = *this; result += change; return result; }
483  inline Self operator- (int change) const { return *this + (-change); }
484 
485  //pointer subtraction
486  inline ptrdiff_t operator- (Self const & other) const { return count_ - other.count_; }
487 
488  //offset dereference
489  inline AtomicType & operator[](int change) { Self result = *this; result += change; return *result; }
490 
491  //comparisons
492  inline bool operator==(Self const & other) const { return current_ == other.current_; }
493  inline bool operator!=(Self const & other) const { return current_ != other.current_; }
494  inline bool operator< (Self const & other) const { return current_ < other.current_; }
495  inline bool operator> (Self const & other) const { return current_ > other.current_; }
496  inline bool operator<=(Self const & other) const { return current_ <= other.current_; }
497  inline bool operator>=(Self const & other) const { return current_ >= other.current_; }
498 
499  IntVec location() const{ return IntVec(xIndex_,yIndex_,zIndex_); }
500 
501  private:
502  AtomicType * current_;
503  int count_;
504  int xIndex_, yIndex_, zIndex_;
505  int yStep_, zStep_;
506  int xExtent_, yExtent_, zExtent_, xyExtent_;
507  };
508 
514  template<typename FieldType>
515  class ConstFieldIterator : public std::iterator<std::random_access_iterator_tag, typename FieldType::value_type> {
517  typedef typename FieldType::value_type AtomicType;
518 
519  public:
521  : current_(NULL),
522  count_(0),
523  xIndex_(0), yIndex_(0), zIndex_(0),
524  yStep_(0), zStep_(0),
525  xExtent_(0), yExtent_(0), zExtent_(0),
526  xyExtent_(0)
527  {}
528 
529  ConstFieldIterator(const AtomicType * const field_values,
530  const MemoryWindow & w)
531  : current_(field_values +
532  w.offset(0) * 1 +
533  w.offset(1) * w.glob_dim(0) +
534  w.offset(2) * w.glob_dim(0) * w.glob_dim(1)),
535  count_(0),
536  xIndex_(0), yIndex_(0), zIndex_(0),
537  yStep_(w.glob_dim(0) - w.extent(0)),
538  zStep_((w.glob_dim(1) - w.extent(1)) * w.glob_dim(0)),
539  xExtent_(w.extent(0)),
540  yExtent_(w.extent(1)),
541  zExtent_(w.extent(2)),
542  xyExtent_(w.extent(0) * w.extent(1))
543  {}
544 
546  : current_(it.current_),
547  count_(it.count_),
548  xIndex_(it.xIndex_),
549  yIndex_(it.yIndex_),
550  zIndex_(it.zIndex_),
551  yStep_(it.yStep_),
552  zStep_(it.zStep_),
553  xExtent_(it.xExtent_),
554  yExtent_(it.yExtent_),
555  zExtent_(it.zExtent_),
556  xyExtent_(it.xyExtent_)
557  {}
558 
559  //immutable dereference
560  inline AtomicType const & operator*() const
561  {
562 # ifndef NDEBUG
563  if( count_ != (xIndex_ +
564  yIndex_ * xExtent_ +
565  zIndex_ * xyExtent_) ){
566  std::ostringstream msg;
567  msg << __FILE__ << " : " << __LINE__ << std::endl
568  << "iterator's internal count is off";
569  throw std::runtime_error(msg.str());
570  }
571  if( xIndex_ >= xExtent_ ||
572  yIndex_ >= yExtent_ ||
573  zIndex_ >= zExtent_ ||
574  xIndex_ < 0 ||
575  yIndex_ < 0 ||
576  zIndex_ < 0 ){
577  std::ostringstream msg;
578  msg << __FILE__ << " : " << __LINE__ << std::endl
579  << "iterator is in an invalid state for dereference";
580  throw std::runtime_error(msg.str());
581  }
582 # endif
583  return *current_;
584  }
585 
586  //increment
587  inline Self & operator++() {
588  current_++; //xStep
589  count_++;
590  xIndex_++;
591  if( xIndex_ == xExtent_ ){
592  current_ += yStep_; //yStep
593  xIndex_ = 0;
594  yIndex_++;
595  if( yIndex_ == yExtent_ ){
596  current_ += zStep_; //zStep
597  yIndex_ = 0;
598  zIndex_++;
599  }
600  }
601  return *this;
602  }
603  inline Self operator++(int) { Self result = *this; ++(*this); return result; }
604 
605  //decrement
606  inline Self & operator--() {
607  --current_; //xStep
608  --count_;
609  --xIndex_;
610  if( xIndex_ == -1 ){
611  current_ -= yStep_; //yStep
612  xIndex_ = xExtent_ - 1;
613  yIndex_--;
614  if( yIndex_ == -1 ){
615  current_ -= zStep_; //zStep
616  yIndex_ = yExtent_ - 1;
617  --zIndex_;
618  }
619  }
620  return *this;
621  }
622  inline Self operator--(int) { Self result = *this; --(*this); return result; }
623 
624  //compound assignment
625  inline Self & operator+=(int change) {
626  //small change (only changes xIndex_)
627  if( (change == 0) || //no change
628  (change > 0 && //positive change
629  change < xExtent_ - xIndex_) ||
630  (change < 0 && //negative change
631  - change < xIndex_) ){
632  current_ += change;
633  xIndex_ += change;
634  count_ += change;
635  }
636  //bigger change (changes yIndex_ and/or zIndex_)
637  else {
638  int new_count = count_ + change;
639  int old_count = count_;
640  current_ += (change + //xStep
641  yStep_ * ((new_count / xExtent_) - (old_count / xExtent_)) +
642  zStep_ * ((new_count / xyExtent_) - (old_count /xyExtent_)));
643  count_ += change;
644  xIndex_ = count_ % xExtent_;
645  yIndex_ = (count_ % xyExtent_) / xExtent_;
646  zIndex_ = count_ / xyExtent_;
647  }
648  return *this;
649  }
650  inline Self & operator-=(int change) { return *this += -change; }
651 
652  //addition/subtraction
653  inline Self operator+ (int change) const { Self result = *this; result += change; return result; }
654  inline Self operator- (int change) const { return *this + (-change); }
655 
656  //iterator subtraction
657  inline ptrdiff_t operator- (Self const & other) const { return count_ - other.count_; }
658 
659  //offset dereference
660  inline AtomicType & operator[](int change) { Self result = *this; result += change; return *result; }
661 
662  IntVec location() const{ return IntVec(xIndex_,yIndex_,zIndex_); }
663 
664  //comparisons
665  inline bool operator==(Self const & other) const { return current_ == other.current_; }
666  inline bool operator!=(Self const & other) const { return current_ != other.current_; }
667  inline bool operator< (Self const & other) const { return current_ < other.current_; }
668  inline bool operator> (Self const & other) const { return current_ > other.current_; }
669  inline bool operator<=(Self const & other) const { return current_ <= other.current_; }
670  inline bool operator>=(Self const & other) const { return current_ >= other.current_; }
671 
672  private:
673  const AtomicType * current_;
674  int count_;
675  int xIndex_, yIndex_, zIndex_;
676  int yStep_, zStep_;
677  int xExtent_, yExtent_, zExtent_, xyExtent_;
678  };
679 
680 } // namespace SpatialOps
681 
682 #endif // SpatialOps_MemoryWindow_h
MemoryWindow remove_ghosts(GhostData const &ghosts) const
return the current memory window with given ghosts removed
Definition: MemoryWindow.h:214
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
MemoryWindow get_window_with_ghost(const IntVec &localDim, const GhostData &ghost, const BoundaryCellInfo &bc)
Obtain the memory window for a field on a patch that is a single, contiguous memory block...
Definition: MemoryWindow.h:314
MemoryWindow(const int npts[3], const int offset[3], const int extent[3])
construct a MemoryWindow object
int flat_index(IntVec loc) const
given the local ijk location (0-based on the local window), obtain the flat index in the global memor...
Definition: MemoryWindow.h:121
std::string print() const
Writes the internals of MemoryWindow to a string.
Definition: MemoryWindow.h:258
bool sanity_check() const
performs basic sanity checks to see if there is anything obviously wrong with this window...
size_t glob_npts() const
obtain the global number of points in the field. Note that this is not necessarily contiguous memory ...
Definition: MemoryWindow.h:155
MemoryWindow reset_ghosts(GhostData const &oldGhosts, GhostData const &newGhosts) const
return the current memory window with oldGhosts removed and newGhosts added
Definition: MemoryWindow.h:198
bool fits_between(MemoryWindow const &inner, MemoryWindow const &outer) const
return if current window fits between the two given windows
Definition: MemoryWindow.h:247
int has_extra(const int dir) const
obtain the number of extra cells actually present on this field due to presence of physical boundarie...
size_t local_npts() const
obtain the local number of points in the field. Note that this is not necessarily contiguous memory ...
Definition: MemoryWindow.h:161
IntVec get_plus() const
obtain the IntVec containing the number of ghost cells on the (+) faces
Definition: GhostData.h:145
int get_dim_with_ghost(const int nNoGhost, const int minusGhost, const int plusGhost, const int bc)
obtain the number of points in the x direction
Definition: MemoryWindow.h:290
IntVec ijk_index_from_global(const int loc) const
given the local flat location (0-based on the global field), obtain the ijk index in the global memor...
Definition: MemoryWindow.h:138
Holds information about the number of ghost cells on each side of the domain.
Definition: GhostData.h:54
IntVec ijk_index_from_local(const int loc) const
given the local flat location (0-based on the local window), obtain the ijk index in the global memor...
Definition: MemoryWindow.h:147
bool operator==(const MemoryWindow &w) const
compare two MemoryWindows for equality
Definition: MemoryWindow.h:180
provides iterator support for SpatialField. Only works for CPU.
Definition: MemoryWindow.h:334
bool fits_in(MemoryWindow const &outer) const
return if the current window fits in given window
Definition: MemoryWindow.h:227
provides iterator support for SpatialField. Only works for CPU.
Definition: MemoryWindow.h:326
IntVec flat_to_ijk(const IntVec &dim, const int pt)
Definition: IntVec.h:76
IntVec get_minus() const
obtain the IntVec containing the number of ghost cells on the (-) faces
Definition: GhostData.h:135
Provides information about boundary cells for various fields.
int ijk_to_flat(const IntVec &dim, const IntVec &loc)
Definition: IntVec.h:62