Array.hpp

/*------------------------------------------------------------------------*/
/*      phdMesh : Parallel Heterogneous Dynamic unstructured Mesh         */
/*                Copyright (2008) Sandia Corporation                     */
/*                                                                        */
/*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
/*  license for use of this work by or on behalf of the U.S. Government.  */
/*                                                                        */
/*  This library is free software; you can redistribute it and/or modify  */
/*  it under the terms of the GNU Lesser General Public License as        */
/*  published by the Free Software Foundation; either version 2.1 of the  */
/*  License, or (at your option) any later version.                       */
/*                                                                        */
/*  This library is distributed in the hope that it will be useful,       */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of        */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU     */
/*  Lesser General Public License for more details.                       */
/*                                                                        */
/*  You should have received a copy of the GNU Lesser General Public      */
/*  License along with this library; if not, write to the Free Software   */
/*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307   */
/*  USA                                                                   */
/*------------------------------------------------------------------------*/

#ifndef util_Array_hpp
#define util_Array_hpp

#ifdef ARRAY_BOUNDS_CHECKING
#define ARRAY_CHECK( X ) X
#else
#define ARRAY_CHECK( X )
#endif

#include <vector>
#include <string>
#include <util/SimpleArrayOps.hpp>

namespace phdmesh {

//----------------------------------------------------------------------
enum ArrayOrder {
  NaturalOrder ,

  FortranOrder ,

  RankZero
};

template< typename Scalar , ArrayOrder Order , 
          class Tag1 = void , class Tag2 = void ,
          class Tag3 = void , class Tag4 = void ,
          class Tag5 = void , class Tag6 = void ,
          class Tag7 = void , class Tag8 = void >
class Array ;

//----------------------------------------------------------------------

template< class ArrayType , class Tag > struct ArrayAppend ;

//----------------------------------------------------------------------
struct ArrayDimTag {

  virtual const char * name() const = 0 ;

  virtual std::string to_string( unsigned dimension ,
                                 unsigned index ) const ;

  virtual unsigned to_index( unsigned dimension ,
                             const std::string & label ) const ; 
 
protected:
  virtual ~ArrayDimTag();
  ArrayDimTag() {}
private:
  ArrayDimTag( const ArrayDimTag & );
  ArrayDimTag & operator = ( const ArrayDimTag & );
};

struct ArrayDimension : public ArrayDimTag {

  const char * name() const ;

  static const ArrayDimension & tag(); 

private:
  ~ArrayDimension();
  ArrayDimension() {}
  ArrayDimension( const ArrayDimension & );
  ArrayDimension & operator = ( const ArrayDimension & );
};

} // namespace phdmesh

//----------------------------------------------------------------------
//----------------------------------------------------------------------

#include <util/ArrayPrivate.hpp>

namespace phdmesh {

//----------------------------------------------------------------------
template< typename Scalar , ArrayOrder array_order ,
          class Tag1 , class Tag2 , class Tag3 , class Tag4 ,
          class Tag5 , class Tag6 , class Tag7 , class Tag8 >
class Array
{
private:
  typedef
    Array<void,array_order,Tag1,Tag2,Tag3,Tag4,Tag5,Tag6,Tag7,Tag8>
      md_type ;
public:
  typedef Scalar  value_type ;

  typedef unsigned size_type ;

  typedef const ArrayDimTag * tag_type ;

  //----------------------------------

  enum { Rank = md_type::Rank };

  enum { Natural = NaturalOrder == array_order };

  enum { Reverse = FortranOrder == array_order };

  enum { Contiguous = true };

  unsigned rank()   const { return Rank ; }

  bool natural()    const { return Natural ; }

  bool reverse()    const { return Reverse ; }

  bool contiguous() const { return Contiguous ; }

  //----------------------------------

#ifndef DOXYGEN_COMPILE

  template < unsigned ordinate >
  struct Tag { typedef typename ArrayTagAt<Array,ordinate>::type type ; };
#endif

  tag_type tag( const unsigned ordinate ) const
    {
      array_check_ordinal( Rank , ordinate );
      return
        array_dim_tags<Tag1,Tag2,Tag3,Tag4,Tag5,Tag6,Tag7,Tag8>()[ordinate];
    }

  //----------------------------------
  template < unsigned ordinate > unsigned dimension() const
    {
      array_check_ordinal_is_less<ordinate,Rank>();
      return ArrayStrideDim<array_order,Rank,ordinate>::dimension(m_stride);
    }

  unsigned dimension( const unsigned ordinate ) const
    {
      array_check_ordinal( Rank , ordinate );
      return ArrayStrideDim<array_order,Rank>::dimension(m_stride,ordinate);
    }

  void dimensions( std::vector<unsigned> & n )
    {
      n.resize( Rank );
      for ( unsigned i = 0 ; i < Rank ; ++i ) { n[i] = dimension(i); }
    }

  size_type size() const { return m_stride[ Rank - 1 ]; }

  //----------------------------------
  typedef typename ArrayTruncate<Array>::type TruncateType ;

  TruncateType truncate( const unsigned i ) const
    {
      TruncateType tmp ;
      tmp.m_ptr = m_ptr + i * ( 1 < Rank ? m_stride[ Rank - 2 ] : 1 );
      Copy<Rank-1>( tmp.m_stride , m_stride );
      return tmp ;
    }

  //----------------------------------
  value_type * contiguous_data() const { return m_ptr ; }

  value_type & operator[]( size_type i ) const
    {
      ARRAY_CHECK( array_check_offset(size(),i) );
      return m_ptr[ i ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ,
                           const unsigned i7 , const unsigned i8 ) const
    { return m_ptr[
        array_offset<array_order,Rank>(m_stride,i1,i2,i3,i4,i5,i6,i7,i8) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ,
                           const unsigned i7 ) const
    { return m_ptr[
        array_offset<array_order,Rank>(m_stride,i1,i2,i3,i4,i5,i6,i7) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ) const
    { return m_ptr[
        array_offset<array_order,Rank>(m_stride,i1,i2,i3,i4,i5,i6) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 ) const
    { return m_ptr[ array_offset<array_order,Rank>(m_stride,i1,i2,i3,i4,i5) ]; }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ) const
    { return m_ptr[ array_offset<array_order,Rank>(m_stride,i1,i2,i3,i4) ]; }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 ) const
    { return m_ptr[ array_offset<array_order,Rank>(m_stride,i1,i2,i3) ]; }

  value_type & operator()( const unsigned i1 , const unsigned i2 ) const
    { return m_ptr[ array_offset<array_order,Rank>(m_stride,i1,i2) ]; }

  value_type & operator()( const unsigned i1 ) const
    { return m_ptr[ array_offset<array_order,Rank>(m_stride,i1) ]; }

  //----------------------------------
  typedef typename ArrayReverse< Array >::type ReverseType ;

  Array() : m_ptr(NULL) { Copy<Rank>( m_stride , (size_type) 0 ); }

  Array( const Array & rhs )
    : m_ptr( rhs.m_ptr ) { Copy<Rank>( m_stride , rhs.m_stride ); }

  Array & operator = ( const Array & rhs )
    {
      m_ptr = rhs.m_ptr ;
      Copy<Rank>( m_stride , rhs.m_stride );
      return *this ;
    }

  Array( const ReverseType & rhs )
    : m_ptr( rhs.m_ptr ) { Copy<Rank>( m_stride , rhs.m_stride ); }

  Array & operator = ( const ReverseType & rhs )
    {
      m_ptr = rhs.m_ptr ;
      Copy<Rank>( m_stride , rhs.m_stride );
      return *this ;
    }

  Array( value_type * arg_ptr , const unsigned * const dims )
    : m_ptr( arg_ptr ) { md_type::assign( m_stride , dims ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 , const unsigned n4 ,
         const unsigned n5 , const unsigned n6 ,
         const unsigned n7 , const unsigned n8 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 , n4 , n5 , n6 , n7 , n8 ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 , const unsigned n4 ,
         const unsigned n5 , const unsigned n6 ,
         const unsigned n7 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 , n4 , n5 , n6 , n7 ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 , const unsigned n4 ,
         const unsigned n5 , const unsigned n6 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 , n4 , n5 , n6 ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 , const unsigned n4 ,
         const unsigned n5 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 , n4 , n5 ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 , const unsigned n4 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 , n4 ); }

  Array( value_type * arg_ptr ,
         const unsigned n1 , const unsigned n2 ,
         const unsigned n3 )
    : m_ptr( arg_ptr )
    { md_type::assign( m_stride , n1 , n2 , n3 ); }

  Array( value_type * arg_ptr , const unsigned n1 , const unsigned n2 )
    : m_ptr( arg_ptr ) { md_type::assign( m_stride , n1 , n2 ); }

  Array( value_type * arg_ptr , const unsigned n1 )
    : m_ptr( arg_ptr ) { md_type::assign( m_stride , n1 ); }

  Array( value_type * arg_ptr )
    : m_ptr( arg_ptr ) { md_type::assign( m_stride ); }

protected:

  value_type * m_ptr ;

  size_type m_stride[ Rank ];

  template< typename , ArrayOrder ,
            class , class , class , class ,
            class , class , class , class >
  friend class phdmesh::Array ;
};

//----------------------------------------------------------------------

#ifndef DOXYGEN_COMPILE

template< typename Scalar >
class Array<Scalar,RankZero,void,void,void,void,void,void,void,void>
{
public:

  typedef Scalar              value_type ;
  typedef unsigned            size_type ;
  typedef const ArrayDimTag * tag_type ;

  //----------------------------------

  enum { Rank       = 0 };
  enum { Natural    = false };
  enum { Reverse    = false };
  enum { Contiguous = true };

  unsigned rank()   const { return Rank ; }
  bool natural()    const { return Natural ; }
  bool reverse()    const { return Reverse ; }
  bool contiguous() const { return Contiguous ; }

  //----------------------------------

  size_type size() const { return m_ptr ? 1 : 0 ; }

  //----------------------------------
  value_type * contiguous_data() const { return m_ptr ; }

  value_type & operator()() const { return *m_ptr ; }

  //----------------------------------
  // Required constructors and assignment operators:

  Array() : m_ptr(NULL) {}

  Array( const Array & rhs ) : m_ptr( rhs.m_ptr ) {}

  Array & operator = ( const Array & rhs )
    { m_ptr = rhs.m_ptr ; return *this ; }

  //----------------------------------
  // Class specific constructors:

  Array( value_type * arg_ptr ) : m_ptr( arg_ptr ) {}

protected:

  value_type * m_ptr ;

  template< typename , ArrayOrder ,
            class , class , class , class ,
            class , class , class , class >
  friend class phdmesh::Array ;
};

#endif /* DOXYGEN_COMPILE */
//----------------------------------------------------------------------
//----------------------------------------------------------------------
template< typename Scalar , ArrayOrder array_order >
class Array<Scalar,array_order,void,void,void,void,void,void,void,void>
{
public:
  typedef Scalar  value_type ;

  typedef unsigned size_type ;

  typedef const ArrayDimTag * tag_type ;

  //----------------------------------

  enum { Natural = NaturalOrder == array_order };

  enum { Reverse = FortranOrder == array_order };

  enum { Contiguous = true };

  unsigned rank()   const { return m_rank ; }

  bool natural()    const { return Natural ; }

  bool reverse()    const { return Reverse ; }

  bool contiguous() const { return Contiguous ; }

  //----------------------------------

  tag_type tag( const unsigned ordinal ) const
    {
      array_check_ordinal( m_rank , ordinal );
      const int i = Natural ? ( m_rank - 1 ) - ordinal : ordinal ;
      return m_tag[i];
    }

  //----------------------------------

  unsigned dimension( const unsigned ordinal ) const
    {
      array_check_ordinal( m_rank , ordinal );
      const int i = Natural ? ( m_rank - 1 ) - ordinal : ordinal ;
      return i ? m_stride[i] / m_stride[i-1] : m_stride[i] ;
    }

  void dimensions( std::vector<unsigned> & n )
    {
      n.resize( m_rank );
      for ( unsigned i = 0 ; i < m_rank ; ++i ) { n[i] = dimension(i); }
    }

  size_type size() const { return m_stride[ m_rank - 1 ]; }

  //----------------------------------
  Array truncate( const unsigned i ) const
    {
      Array tmp ;
      if ( 1 < m_rank ) {
        tmp.m_ptr  = m_ptr + m_stride[ m_rank - 2 ] * i ;
        tmp.m_rank = m_rank - 1 ;
        unsigned k ;
        for ( k = 0 ; k < m_rank - 1 ; ++k ) { tmp.m_stride[i] = m_stride[i] ; }
        for (       ; k < 8          ; ++k ) { tmp.m_stride[i] = 0 ; }
        for ( k = 0 ; k < m_rank - 1 ; ++k ) { tmp.m_tag[i] = m_tag[i] ; }
        for (       ; k < 8          ; ++k ) { tmp.m_tag[i] = NULL ; }
      }
      return tmp ;
    }

  value_type * contiguous_data() const { return m_ptr ; }

  value_type & operator[]( size_type i ) const
    {
      ARRAY_CHECK( array_check_offset(size(),i) );
      return m_ptr[ i ];
    }

  //----------------------------------
  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ,
                           const unsigned i7 , const unsigned i8 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 8 ) );
      return m_ptr[
        array_offset<array_order,8>(m_stride,i1,i2,i3,i4,i5,i6,i7,i8) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ,
                           const unsigned i7 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 7 ) );
      return m_ptr[
        array_offset<array_order,7>(m_stride,i1,i2,i3,i4,i5,i6,i7) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 , const unsigned i6 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 6 ) );
      return m_ptr[ array_offset<array_order,6>(m_stride,i1,i2,i3,i4,i5,i6) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ,
                           const unsigned i5 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 5 ) );
      return m_ptr[ array_offset<array_order,5>(m_stride,i1,i2,i3,i4,i5) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 , const unsigned i4 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 4 ) );
      return m_ptr[ array_offset<array_order,4>(m_stride,i1,i2,i3,i4) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ,
                           const unsigned i3 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 3 ) );
      return m_ptr[ array_offset<array_order,3>(m_stride,i1,i2,i3) ];
    }

  value_type & operator()( const unsigned i1 , const unsigned i2 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 2 ) );
      return m_ptr[ array_offset<array_order,2>(m_stride,i1,i2) ];
    }

  value_type & operator()( const unsigned i1 ) const
    {
      ARRAY_CHECK( array_check_rank( m_rank , 1 ) );
      return m_ptr[ array_offset<array_order,1>(m_stride,i1) ];
    }

  //----------------------------------
  typedef typename ArrayReverse< Array >::type ReverseType ;

  //----------------------------------

  Array()
    : m_ptr(NULL), m_rank(0)
    {
      Copy<8>( m_stride , (size_type) 0 );
      Copy<8>( m_tag , (tag_type) NULL );
    }

  Array( const Array & rhs )
    : m_ptr( rhs.m_ptr ), m_rank( rhs.m_rank )
    {
      Copy<8>( m_stride , rhs.m_stride );
      Copy<8>( m_tag , rhs.m_tag );
    }

  Array & operator = ( const Array & rhs )
    {
      m_ptr = rhs.m_ptr ;
      m_rank = rhs.m_rank ;
      Copy<8>( m_stride , rhs.m_stride );
      Copy<8>( m_tag , rhs.m_tag );
      return *this ;
    }

  Array( const ReverseType & rhs )
    : m_ptr( rhs.m_ptr ), m_rank( rhs.m_rank )
    {
      Copy<8>( m_stride , rhs.m_stride );
      Copy<8>( m_tag , rhs.m_tag );
    }

  Array & operator = ( const ReverseType & rhs )
    {
      m_ptr = rhs.m_ptr ;
      m_rank = rhs.m_rank ;
      Copy<8>( m_stride , rhs.m_stride );
      Copy<8>( m_tag , rhs.m_tag );
      return *this ;
    }

  //----------------------------------

  template< ArrayOrder order ,
            class Tag1 , class Tag2 , class Tag3 , class Tag4 ,
            class Tag5 , class Tag6 , class Tag7 , class Tag8 >
  Array(
    const Array<value_type,order,Tag1,Tag2,Tag3,Tag4,Tag5,Tag6,Tag7,Tag8> & rhs )
  : m_ptr( rhs.m_ptr ), m_rank( 0 )
  {
    typedef Array<value_type,order,Tag1,Tag2,Tag3,Tag4,Tag5,Tag6,Tag7,Tag8> a_t ;
    enum { inRank    = a_t::Rank };
    enum { inNatural = a_t::Natural };
    m_rank = inRank ;
    Copy< inRank >(     m_stride , rhs.m_stride );
    Copy< 8 - inRank >( m_stride + inRank , (size_type) 0 );
    unsigned i = 0 ;
    if ( inNatural ) {
      for ( ; i < inRank ; ++i ) { m_tag[i] = rhs.tag((inRank-1)-i); }
    }
    else {
      for ( ; i < inRank ; ++i ) { m_tag[i] = rhs.tag(i); }
    }
    for ( ; i < 8 ; ++i ) { m_tag[i] = NULL ; }
  }

  //----------------------------------
  // Class specific constructors:

  Array( value_type * ptr ,
         const unsigned rank ,
         const unsigned * const dims ,
         const tag_type  * const tags )
    : m_ptr( ptr ), m_rank( rank )
    {
      if ( Natural ) {
        size_type n = 1 ;
        unsigned i ;
        for ( i = 0 ; i < rank ; ++i ) { m_stride[i] = n *= dims[(rank-1)-i]; }
        for (       ; i < 8    ; ++i ) { m_stride[i] = 0 ; }
        for ( i = 0 ; i < rank ; ++i ) { m_tag[i] = tags[(rank-1)-i]; }
        for (       ; i < 8    ; ++i ) { m_tag[i] = NULL ; }
      }
      else {
        size_type n = 1 ;
        unsigned i ;
        for ( i = 0 ; i < rank ; ++i ) { m_stride[i] = n *= dims[i] ; }
        for (       ; i < 8    ; ++i ) { m_stride[i] = 0 ; }
        for ( i = 0 ; i < rank ; ++i ) { m_tag[i] = tags[i]; }
        for (       ; i < 8    ; ++i ) { m_tag[i] = NULL ; }
      }
    }

protected:

  value_type * m_ptr ;

  unsigned     m_rank ;

  size_type    m_stride[8];

  tag_type     m_tag[8] ;

  template< typename , ArrayOrder ,
            class , class , class , class ,
            class , class , class , class >
  friend class phdmesh::Array ;
};

//----------------------------------------------------------------------
//----------------------------------------------------------------------

} // namespace phdmesh

#undef ARRAY_CHECK

#endif