#ifndef __GRID_H__
#define __GRID_H__

#include <string>
#include <vector>

#include "vtkSmartPointer.h"

class vtkUnstructuredGrid;
class vtkCPDataDescription;
class vtkDataObject;

namespace CamAdaptor
{

enum GridType
{
  RECTILINEAR,
  SPHERE,
  CUBE_SPHERE,
  GRID_TYPE_COUNT
};



/// Creates and accumulates data for a 2D and a 3D grid. It generates either
/// a parallelepiped or sphere for the Finite Volume (FV) dynamic core
/// or a cube sphere for the Structured Element (SE) dynamic core.
/// Grids are stored as unstructured grids, and are reconstructed from points.
template<GridType gridType>
class Grid
{
public:

  Grid();
  /// Deletes data used to build the grids. Note that the grid memory
  /// is managed by the Catalyst Coprocessor.
  ~Grid();

  /// Creates a 2D and a 3D grid
  void Create();
  void SetMpiRank(int rank)
  {
    this->MpiRank = rank;
  }
  void SetChunkCapacity(int capacity)
  {
    this->ChunkCapacity = capacity;
  }
  void SetNCells2d(int ncells)
  {
    this->NCells2d = ncells;
  }
  void SetNLon(int nlon)
  {
    this->NLon = nlon;
  }
  void SetNLat(int nlat)
  {
    this->NLat = nlat;
  }
  void SetLev(int nlev, double* lev)
  {
    this->NLev = nlev;
    if (lev)
      {
      this->Lev = new double[nlev];
      std::copy(lev, lev + nlev, this->Lev);
      }
  }
  /// Adds the points and the cells for a vertical column to the grid 
  void AddPointsAndCells(double lonRad, double latRad);

  /// Sets attributes for a chunk (a list of vertical columns) to 
  /// the 2D and 3D grids
  void SetAttributeValue(int chunkSize,
                         double* lonRad, double* latRad,
                         double* psScalar, double *tScalar, double* uScalar, 
                         double* vScalar);

  /// Returns the 2D grid
  vtkSmartPointer<vtkUnstructuredGrid> GetGrid2d() const
  {
    return this->Grid2d;
  }
  /// Returns the 3D grid
  vtkSmartPointer<vtkUnstructuredGrid> GetGrid3d() const
  {
    return this->Grid3d;
  }
  /// Attach the grid to the coprocessor data
  static bool SetToCoprocessor(
    vtkCPDataDescription* coprocessorData, const char* name,
    vtkSmartPointer<vtkUnstructuredGrid> grid);


  ///////////////////////////////////////////////////////////////////////////
  /// used for the FV dynamic core only
  void SetLonStep(int step)
  {
    this->LonStep = step;
  }
  void SetLatStep(int step)
  {
    this->LatStep = step;
  }
  /// used for the SE dynamic core only
  void SetCubeGridPoints(
    int ne, int np, int lonSize, double* lonRad, int latSize, double* latRad);
  ///////////////////////////////////////////////////////////////////////////

private:
  /// Print routine for debugging
  void PrintAddChunk(
    int* nstep, int* chunkCols,
    double* lonRad, double* latRad, double* psScalar, double *tScalar);

  /// Creates the 2D grid and the data used to add attributes to the grid
  vtkSmartPointer<vtkUnstructuredGrid> CreateGrid2d();

  /// Creates the 3D grid and the data used to add attributes to the grid
  vtkSmartPointer<vtkUnstructuredGrid> CreateGrid3d();

  /// The level for a grid measures preasure which decreases with height.
  /// This transformation sets 
  /// - max pressure to 0 (surface level) and min pressure to 1 (SPHERE)
  /// - negates the pressure so that globe surface (high pressure) is closer to 
  ///   origin than high in the atmosphere (low pressure) (RECTILINEAR)
  /// These transformations match the transformation in the history file.
  double LevelToRadius(double level) const;

  /// Compute the  location of  the points surounding  a cell  at index
  /// 'levIndex', 'latIndex' or 'lonIndex'
  void GetLevMinusPlus(int levIndex, double* levMinus, double* levPlus) const;
  void GetLatMinusPlus(int latIndex, double* latMinus, double* latPlus) const;
  void GetLonMinusPlus(int lonIndex, double* lonMinus, double* lonPlus) const;

  /// Returns (lonIndex, latIndex) for FV dynamic core (regular grid) or 
  /// (s, t, face) for SE dynamic core (cube sphere grid).
  void GetValueIndex(double lonRad, double latRad,
                     std::vector<std::vector<int> >& index) const;

  /// Points for a 2D cell surrounding value at index for 2D or (index, level) for 3D.
  void GetCellPoints(const int index[3], double cellPoints[4][3]) const;
  void GetCellPoints(const int index[3], int level, double cell[8][3]) const;

  /// Returns unique ids for cells' points associated with value at 'index'
  void GetPointIds(const int index[3], int pointIndexes[4]) const;
  void GetPointIds(const int index[3], int level, int pointIndexes[8]) const;

  /// Returns unique id for a cell at index
  int GetCellId(const int index[3]) const;
  int GetCellId(const int index[3], int level) const;

  /// converts from a (x,y,z) coordinate on the cube to an (s, t, face) index
  /// index is between (0,0,0) - (NLon-1, NLon-1, NLon-1). Returns the
  /// index in the vector that has the minimum face.
  size_t CubeToIndex (double side, const std::vector<double>& cubeCoordinates, double v[3],
                      std::vector<std::vector<int> >& index) const;

  /// converts from (s, t, face) to cartesian coordinates on the cube
  /// the index is double[3] because we can get indexes in between integer positions.
  void IndexToCube (const std::vector<double>& cubeCoordinates,
                    double index[3], double v[3]) const;
  void CubeToSpherical (double v[3], double* lonRad, double* latRad) const;

  /// converts 'index' to an index where the face is minimum. A conversion
  /// is made only for indexes that are shared between 2 or 3 faces.
  void MinFaceIndex(int side, int index[3]) const;

private:
  int MpiRank;        /// MPI rank
  int ChunkCapacity;  /// maximum number of (vertical) columns in a chunk
  int NCells2d;       /// total number of 2D cells on a MPI processor
  int NLon, NLat, NLev; /// number of unique longitude, latitude or level values.
                        /// for SE, NLon == NLat == number of values on the side
                        /// of the cube (number of elements + 1)
  double* Lev;        /// level values
  ///////////////////////////////////////////////////////////////////////////
  /// used for the SE dynamic core only
  double Radius;      /// radius around a cube with side ne * np.
  int Ne;             /// the number of spectral  elements ne along the
                      /// edge of each cube face
  int Np;             /// each element is further subdivided by the Np x Np
                      /// Gauss-Lobatto-Legendre (GLL) quadrature
                      /// points where GLL stands for
  /// used for FV dynamic core only
  double LonStep;     /// longitude step in degrees
  double LatStep;     /// latitude step in degrees
  ///////////////////////////////////////////////////////////////////////////
  class Internals;
  Internals* Impl;

  /// 2d grid
  vtkSmartPointer<vtkUnstructuredGrid> Grid2d;//the grid
  int RankArrayIndex2d;       //index of the rank array
  int CoordArrayIndex2d;      //index of the coordinates array (lon x lat)
  int PSArrayIndex;           /// index of the PS array

  /// 3d grid
  vtkSmartPointer<vtkUnstructuredGrid> Grid3d;/// the grid
  int RankArrayIndex3d;       /// index of the rank array
  int CoordArrayIndex3d;      /// index of the coordinates array (lon x lat x lev)
  int TArrayIndex;            /// index of the attribute array for T
  int UArrayIndex;            /// index of the attribute array for U
  int VArrayIndex;            /// index of the attribute array for V
};

};
#endif
