GridArray.h

#ifndef UTIL_GRID_ARRAY_H
#define UTIL_GRID_ARRAY_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <util/misc/Memory.h>
#include <util/global.h>

namespace Util
{

   template <typename Data>
   class GridArray
   {

   public:

      GridArray();

      GridArray(const GridArray<Data>& other);

      ~GridArray();

      GridArray<Data>& operator = (const GridArray<Data>& other);

      // Initialization

      void allocate(const IntVector& dimensions);

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

      bool isAllocated() const;

      // Grid Interface

      const IntVector& dimensions();

      int dimension(int i) const;

      int size() const;

      IntVector position(int rank) const;

      int rank(const IntVector& position) const;

      bool isInGrid(int coordinate, int i) const;

      bool isInGrid(IntVector& position) const;

      int shift(int& coordinate, int i) const;

      IntVector shift(IntVector& position) const;

      // Array Interface

      const Data& operator[] (int rank) const;

      Data& operator[] (int rank);

      const Data& operator() (const IntVector& position) const;

      Data& operator() (const IntVector& position);

      /*
      * Return pointer to underlying 1D C-array.
      */
      Data* data();

   private:

      Data*  data_;

      IntVector offsets_;

      IntVector dimensions_;

      int size_;

   };

   // Method definitions

   template <typename Data>
   inline GridArray<Data>::GridArray()
    : data_(0),
      offsets_(),
      dimensions_(),
      size_(0)
   {}

   /*
   * Destructor.
   *
   * Delete dynamically allocated C array.
   */
   template <typename Data>
   GridArray<Data>::~GridArray()
   {
      if (data_) {
         Memory::deallocate<Data>(data_, size_);
         size_ = 0;
      }
   }

   /*
   * Copy constructor.
   */
   template <typename Data>
   GridArray<Data>::GridArray(const GridArray<Data>& other)
    : data_(0),
      offsets_(),
      dimensions_(),
      size_(0)
   {
      // Precondition
      if (other.data_ == 0) {
         UTIL_THROW("Other GridArray must be allocated");
      }
      if (isAllocated()) {
         UTIL_THROW("GridArray already allocated in copy constructor");
      }

      allocate(other.dimensions_);
      if (offsets_ != other.offsets_ ) {
         UTIL_THROW("Unequal offsets");
      }
      if (size_ != other.size_ ) {
         UTIL_THROW("Unequal sizes");
      }
      for (int i = 0; i < size_; ++i) {
         data_[i] = other.data_[i];
      }
   }

   /*
   * Assignment.
   */
   template <typename Data>
   GridArray<Data>& GridArray<Data>::operator = (const GridArray<Data>& other)
   {
      // Check for self assignment.
      if (this == &other) {
         return *this;
      }

      // Precondition
      if (other.data_ == 0) {
         UTIL_THROW("Other GridArray must be allocated before assignment");
      }

      // If this GridArray if not allocated, allocate now.
      // If it is allocated, check that dimensions are equal.
      if (!isAllocated()) {
         allocate(other.dimensions_);
      } else {
         if (dimensions_ != other.dimensions_ ) {
            UTIL_THROW("Unequal dimensions");
         }
      }
      if (offsets_ != other.offsets_ ) {
         UTIL_THROW("Unequal offsets");
      }
      if (size_ != other.size_ ) {
         UTIL_THROW("Unequal sizes");
      }

      // Copy elements
      for (int i = 0; i < size_; ++i) {
         data_[i] = other.data_[i];
      }

      return *this;
   }

   /*
   * Set dimensions and allocate memory.
   */
   template <typename Data>
   void GridArray<Data>::allocate(const IntVector& dimensions)
   {
      if (isAllocated()) {
         UTIL_THROW("Attempt to re-allocate a GridArray");
      }
      for (int i = 0; i < Dimension; ++i) {
         if (dimensions[i] <= 0) {
            UTIL_THROW("Dimension not positive");
         }
      }
      dimensions_ = dimensions;
      offsets_[Dimension -1] = 1;
      for (int i = Dimension - 1; i > 0; --i) {
         offsets_[i-1] = offsets_[i]*dimensions_[i];
      }
      size_ = offsets_[0]*dimensions_[0];
      Memory::allocate<Data>(data_, size_);
   }

   /*
   * Serialize a GridArray to/from an Archive.
   */
   template <class Data>
   template <class Archive>
   void GridArray<Data>::serialize(Archive& ar, const unsigned int version)
   {
      IntVector dimensions;
      if (Archive::is_saving()) {
         dimensions = dimensions_;
      }
      ar & dimensions;
      if (Archive::is_loading()) {
         if (!isAllocated()) {
            allocate(dimensions);
         }
      }
      for (int i = 0; i < size_; ++i) {
         ar & data_[i];
      }
   }

   /*
   * Get IntVector of dimensions.
   */
   template <typename Data>
   inline const IntVector& GridArray<Data>::dimensions()
   {  return dimensions_; }

   /*
   * Get dimension in direction i.
   */
   template <class Data>
   inline int GridArray<Data>::dimension(int i) const
   {  return dimensions_[i]; }

   /*
   * Get total number of grid points.
   */
   template <class Data>
   inline int GridArray<Data>::size() const
   {  return size_; }

   /*
   * Calculate 1D rank from grid coordinates.
   */
   template <typename Data>
   #ifdef UTIL_DEBUG
   int GridArray<Data>::rank(const IntVector& position) const
   {
      int result = 0;
      int i;
      for (i = 0; i < Dimension - 1; ++i) {
         assert(position[i] >= 0);
         assert(position[i] < dimensions_[i]);
         result += position[i]*offsets_[i];
      }
      assert(position[i] >= 0);
      assert(position[i] < dimensions_[i]);
      result += position[i];
      return result;
   }
   #else
   inline int GridArray<Data>::rank(const IntVector& position) const
   {
      return (position[0]*offsets_[0] + position[1]*offsets_[1] + position[2]);
   }
   #endif

   /*
   * Calculate coordinates from 1D rank.
   */
   template <typename Data>
   IntVector GridArray<Data>::position(int rank) const
   {
      IntVector position;
      int remainder = rank;

      int i;
      for (i = 0; i < Dimension - 1; ++i) {
         position[i] = remainder/offsets_[i];
         remainder -= position[i]*offsets_[i];
      }
      position[i] = remainder;
      return position;
   }

   /*
   * Test if a single coordinate is within range of grid.
   */
   template <typename Data>
   bool GridArray<Data>::isInGrid(int coordinate, int i) const
   {
      bool result = true;
      if (coordinate <  0)
         result = false;
      if (coordinate >= dimensions_[i])
         result = false;
      return result;
   }

   /*
   * Test if a IntVector position is within primary grid domain.
   */
   template <typename Data>
   bool GridArray<Data>::isInGrid(IntVector& position) const
   {
      bool result = true;
      for (int i = 0; i < Dimension; ++i) {
         if (position[i] <  0)
            result = false;
         if (position[i] >= dimensions_[i])
            result = false;
      }
      return result;
   }

   /*
   * Shift a 1D coordinate to primary domain.
   */
   template <typename Data>
   int GridArray<Data>::shift(int& coordinate, int i) const
   {
      int shift;
      if (coordinate >= 0) {
         shift = coordinate/dimensions_[i];
      } else {
         shift = -1 + ((coordinate+1)/dimensions_[i]);
      }
      coordinate -= shift*dimensions_[i];
      return shift;
   }

   /*
   * Shift a IntVector position to primary domain.
   */
   template <typename Data>
   IntVector GridArray<Data>::shift(IntVector& position) const
   {
      IntVector shifts;
      for (int i = 0; i < Dimension; ++i) {
         shifts[i] = shift(position[i], i);
      }
      return shifts;
   }

   /*
   * Return element by const reference, indexed by 1D rank.
   */
   template <typename Data>
   inline const Data& GridArray<Data>::operator[] (int rank) const
   {  return *(data_ + rank); }

   /*
   * Return element by reference, indexed by 1D rank.
   */
   template <typename Data>
   inline Data& GridArray<Data>::operator[] (int rank)
   {  return *(data_ + rank); }

   /*
   * Return element by const reference, indexed by IntVector of coordinates
   */
   template <typename Data>
   inline
   const Data& GridArray<Data>::operator() (const IntVector& position) const
   {  return *(data_ + rank(position)); }

   /*
   * Return element by reference, indexed by IntVector of coordinates
   */
   template <typename Data>
   inline Data& GridArray<Data>::operator() (const IntVector& position)
   {  return *(data_ + rank(position)); }

   /*
   * Return pointer to underlying 1D C-array.
   */
   template <typename Data>
   inline Data* GridArray<Data>::data() 
   {  return data_;}

   /*
   * Return true if the GridArray has been allocated, false otherwise.
   */
   template <class Data>
   inline bool GridArray<Data>::isAllocated() const
   {  return (bool)(data_ != 0); }

}
#endif