stencil2D.h

Go to the documentation of this file.

#ifndef _DLRNUMERIC_STENCIL2D_H_
#define _DLRNUMERIC_STENCIL2D_H_

#include <vector>
#include <dlrNumeric/array2D.h>
#include <dlrNumeric/index2D.h>

namespace dlr {

  namespace numeric {
    
    template <class Type, int Size>
    class StencilIterator
    {
    public:

      StencilIterator(Type* ptr, int* incrementPtr, size_t position=0)
        : m_incrementPtr(incrementPtr),
          m_position(position),
          m_ptr(ptr) {}

      bool
      operator!=(const StencilIterator<Type, Size>& other) {
        return m_ptr != other.m_ptr;
      }

    
      bool
      operator==(const StencilIterator<Type, Size>& other) {
        return m_ptr == other.m_ptr;
      }

    
      StencilIterator<Type, Size>&
      operator++() {m_ptr += m_incrementPtr[m_position++]; return *this;}

    
      StencilIterator<Type, Size>
      operator++(int) {
        StencilIterator<Type, Size> result(*this);
        ++(*this);
        return result;
      }


      StencilIterator<Type, Size>&
      operator--() {m_ptr += m_incrementPtr[--m_position]; return *this;}

    
      StencilIterator<Type, Size>
      operator--(int) {
        StencilIterator<Type, Size> result(*this);
        --(*this);
        return result;
      }


      Type&
      operator*() {return *m_ptr;}

    
      Type* operator->() {return m_ptr;}

    private:
      int* m_incrementPtr;
      size_t m_position;
      Type* m_ptr;
    };


    template <class Type, int Size>
    class Stencil2D {
    public:
      Stencil2D();


      Stencil2D(size_t rows, size_t columns);
    

      Stencil2D(const Array2D<bool>& pattern);


      Stencil2D(const std::vector<Index2D>& pattern);
    

      ~Stencil2D() {}


      inline void
      advance() {++m_ptr;}
    

      StencilIterator<Type, Size>
      begin() {return StencilIterator<Type, Size>(
                 m_ptr + m_offsetArray[0], m_incrementArray, 0);}


      StencilIterator<Type, Size>
      end() {return StencilIterator<Type, Size>(
               m_ptr + m_targetSize, m_incrementArray, m_numberOfElements);}
    

      inline Type&
      getReference(size_t elementNumber);


      inline Type
      getValue(size_t elementNumber) const;


      inline void
      goTo(size_t index);


      inline void
      goTo(size_t row, size_t column);


      void
      setPattern(const Array2D<bool>& pattern);


      template <class Type2>
      void
      setTarget(const Array2D<Type2>& target);


      template <class Type2>
      void
      setTarget(Array2D<Type2>& target);
    
      // private:

      inline void
      checkBounds(size_t elementNumber) const;


      // Data members to allow moving around in the target array.
      Type* m_basePtr;
      size_t m_rows;
      size_t m_columns;

      // Data members to allow bounds checking.
      int m_targetSize;

      // Data members which actually do the work.
      size_t m_numberOfElements;
      Type* m_ptr;
      int m_offsetArray[Size];
      int m_incrementArray[Size];
      Index2D m_patternArray[Size];
    };

  } // namespace numeric

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using numeric::StencilIterator;
  using numeric::Stencil2D;

} // namespace dlr


/* ============ Definitions of inline & template functions ============ */


#include <cmath>

namespace dlr {

  namespace numeric {
    
    template <class Type, int Size>
    Stencil2D<Type, Size>::
    Stencil2D()
      : m_basePtr(0),
        m_rows(0),
        m_columns(0),
        m_targetSize(0),
        m_numberOfElements(0),
        m_ptr(0),
        m_offsetArray(),
        m_incrementArray(),
        m_patternArray()
    {
      for(size_t elementIndex = 0; elementIndex < Size; ++elementIndex) {
        m_offsetArray[elementIndex] = 0;
        m_incrementArray[elementIndex] = 0;
        m_patternArray[elementIndex] = Index2D(0, 0);
      }
    }


    template <class Type, int Size>
    Stencil2D<Type, Size>::
    Stencil2D(size_t rows, size_t columns)
      : m_basePtr(0),
        m_rows(0),
        m_columns(0),
        m_targetSize(0),
        m_numberOfElements(rows * columns),
        m_ptr(0),
        m_offsetArray(),
        m_incrementArray(),
        m_patternArray()
    {
      if(rows * columns > Size) {
        std::ostringstream message;
        message << "Requested size, (" << rows << ", " << columns << ") "
                << "is too big for a " << Size << " element stencil.";
        DLR_THROW(ValueException, "Stencil2D::Stencil2D()",
                  message.str().c_str());
      }

      // Record the rows x columns ROI.
      size_t index0 = 0;
      for(int row = 0; row < static_cast<int>(rows); ++row) {
        for(int column = 0; column < static_cast<int>(columns); ++column) {
          m_offsetArray[index0] = 0;
          m_incrementArray[index0] = 0;
          m_patternArray[index0] = Index2D(row, column);
          ++index0;
        }
      }

      // Fill in any unused elements with zeros.
      while(index0 < Size) {
        m_offsetArray[index0] = 0;
        m_incrementArray[index0] = 0;
        m_patternArray[index0] = Index2D(0, 0);
        ++index0;
      }
    }
  
    
    template <class Type, int Size>
    Stencil2D<Type, Size>::
    Stencil2D(const Array2D<bool>& pattern)
      : m_basePtr(0),
        m_rows(0),
        m_columns(0),
        m_targetSize(0),
        m_numberOfElements(0),
        m_ptr(0),
        m_offsetArray(),
        m_incrementArray(),
        m_patternArray()
    {
      this->setPattern(pattern);
    }


    template <class Type, int Size>
    Stencil2D<Type, Size>::
    Stencil2D(const std::vector<Index2D>& pattern)
      : m_basePtr(0),
        m_rows(0),
        m_columns(0),
        m_targetSize(0),
        m_numberOfElements(pattern.size()),
        m_ptr(0),
        m_offsetArray(),
        m_incrementArray(),
        m_patternArray()
    {
      if(pattern.size() > Size) {
        std::ostringstream message;
        message << "Requested pattern has too many elements for a "
                << Size << " element stencil.";
        DLR_THROW(ValueException, "Stencil2D::Stencil2D()",
                  message.str().c_str());
      }
      std::fill(m_offsetArray.begin(), m_offsetArray.end(), 0);
      std::fill(m_incrementArray.begin(), m_incrementArray.end(), 0);
      std::copy(pattern.begin(), pattern.end(), m_patternArray.begin());
      for(size_t index0 = pattern.size(); index0 < Size; ++index0) {
        m_patternArray[index0] = Index2D(0, 0);
      }
    }

  
    template <class Type, int Size>
    Type&
    Stencil2D<Type, Size>::
    getReference(size_t elementNumber)
    {
      this->checkBounds(elementNumber);
      return *(m_ptr + m_offsetArray[elementNumber]);
    }
      

    template <class Type, int Size>
    Type
    Stencil2D<Type, Size>::
    getValue(size_t elementNumber) const
    {
      this->checkBounds(elementNumber);
      return *(m_ptr + m_offsetArray[elementNumber]);
    }


    template <class Type, int Size>
    void
    Stencil2D<Type, Size>::
    goTo(size_t index)
    {
      m_ptr = m_basePtr + index;
    }


    template <class Type, int Size>
    void
    Stencil2D<Type, Size>::
    goTo(size_t row, size_t column)
    {
      m_ptr = m_basePtr + row * m_columns + column;
    }


    template <class Type, int Size>
    void
    Stencil2D<Type, Size>::
    setPattern(const Array2D<bool>& pattern)
    {
      // Record which of the elements in pattern are true.
      size_t index0 = 0;
      for(int row = 0; row < static_cast<int>(pattern.rows()); ++row) {
        for(int column = 0; column < static_cast<int>(pattern.columns());
            ++column) {
          if(pattern(row, column)) {
            if(index0 >= Size) {
              std::ostringstream message;
              message << "Requested pattern has too many elements for a "
                      << Size << " element stencil.";
              DLR_THROW(ValueException, "Stencil2D::Stencil2D()",
                        message.str().c_str());
            }
            m_offsetArray[index0] = 0;
            m_incrementArray[index0] = 0;
            m_patternArray[index0] = Index2D(row, column);
            ++index0;
          }
        }
      }
      m_numberOfElements = index0;

      // Fill in any unused elements with zeros.
      while(index0 < Size) {
        m_offsetArray[index0] = 0;
        m_incrementArray[index0] = 0;
        m_patternArray[index0] = Index2D(0, 0);
        ++index0;
      }
    }

    
    template <class Type, int Size> template <class Type2>
    void
    Stencil2D<Type, Size>::
    setTarget(const Array2D<Type2>& target)
    {
      m_basePtr = target.data();
      m_rows = target.rows();
      m_columns = target.columns();
      m_targetSize = static_cast<int>(target.size());
      m_ptr = m_basePtr;

      for(size_t elementIndex = 0; elementIndex < m_numberOfElements;
          ++elementIndex) {
        Index2D targetIndex = m_patternArray[elementIndex];
        m_offsetArray[elementIndex] =
          static_cast<int>(targetIndex.getColumn() + targetIndex.getRow() * m_columns);
      }
      for(size_t elementIndex = 0; elementIndex < m_numberOfElements - 1;
          ++elementIndex) {
        m_incrementArray[elementIndex] =
          m_offsetArray[elementIndex + 1] - m_offsetArray[elementIndex];
      }
      m_incrementArray[m_numberOfElements - 1] =
        m_targetSize - m_offsetArray[m_numberOfElements - 1];
    }


    template <class Type, int Size> template <class Type2>
    void
    Stencil2D<Type, Size>::
    setTarget(Array2D<Type2>& target)
    {
      m_basePtr = target.data();
      m_rows = target.rows();
      m_columns = target.columns();
      m_targetSize = target.size();
      m_ptr = m_basePtr;

      for(size_t elementIndex = 0; elementIndex < m_numberOfElements;
          ++elementIndex) {
        Index2D targetIndex = m_patternArray[elementIndex];
        m_offsetArray[elementIndex] =
          targetIndex.getColumn() + targetIndex.getRow() * m_columns;
      }
      for(size_t elementIndex = 0; elementIndex < m_numberOfElements - 1;
          ++elementIndex) {
        m_incrementArray[elementIndex] =
          m_offsetArray[elementIndex + 1] - m_offsetArray[elementIndex];
      }
      m_incrementArray[m_numberOfElements - 1] =
        m_targetSize - m_offsetArray[m_numberOfElements - 1];
    }


    template <class Type, int Size>
    inline void
    Stencil2D<Type, Size>::
    checkBounds(size_t elementNumber) const
    {
#ifdef _DLRNUMERIC_CHECKBOUNDS_
      if(elementNumber >= m_numberOfElements) {
        std::ostringstream message;
        message << "Element number " << elementNumber
                << " is invalid for a(n) " << m_numberOfElements
                << " element stencil. ";
        DLR_THROW(IndexException, "Stencil2D::checkBounds(size_t)",
                  message.str().c_str());
      }

      Type* resultPtr = m_ptr + this->m_offsetArray[elementNumber];
      if(resultPtr < m_basePtr) {
        DLR_THROW(IndexException, "Stencil2D::checkBounds(size_t)",
                  "Pointer has retreated too far.  "
                  "Stencil no longer addresses valid memory.");
      }

      if(resultPtr >= m_basePtr + m_targetSize) {
        DLR_THROW(IndexException, "Stencil2D::checkBounds(size_t)",
                  "Pointer has advanced too far.  "
                  "Stencil no longer addresses valid memory.");
      }
#endif 
    }
  
  } // namespace numeric

} // namespace dlr

#endif /* #ifndef _DLRNUMERIC_STENCIL2D_H_ */