#ifndef _GRID_H_
#define _GRID_H_

#include <assert.h>
#include <stdio.h>
#include <GL/gl.h>

#include "argparser.h"
#include "vectors.h"
#include "cell.h"

class ArgParser;

// ========================================================================
// ========================================================================

class Grid {

public:

  // ========================
  // CONSTRUCTOR & DESTRUCTOR
  Grid(ArgParser *_args);
  ~Grid();
  void Load();
  void Initialize();

  // ===========================
  // GENERAL GRID CELL ACCESSORS
  int Index(int x, int y, int z) const {
    assert (x >= 0 && x < nx);
    assert (y >= 0 && y < ny);
    assert (z >= 0 && z < nz);
    return x*ny*nz + y*nz + z;
  }

  Cell* getCell(int x, int y, int z) const {
    assert (x >= 0 && x <= nx);
    assert (y >= 0 && y <= ny);
    assert (z >= 0 && z <= nz);
    return &cells[Index(x,y,z)];
  }

  // *********************************************************************  
  // ASSIGNMENT:
  //
  // Clamping the cell indices may or may not be appropriate in different 
  // situations (this could be the source of some bugs).
  //
  // If you want to simulate an "infinite" ocean of water by wrapping 
  // around the x and z dimensions you'll need to change how this is used.
  //
  Cell* getClampedCell(int x, int y, int z) const {
    // clamp at boundaries
    if (x < 0) x = 0;
    if (x >= nx) x = nx-1;
    if (y < 0) y = 0;
    if (y >= ny) y = ny-1;
    if (z < 0) z = 0;
    if (z >= nz) z = nz-1;
    return getCell(x,y,z);
  }
  //
  // *********************************************************************  

  float getPressure(int i, int j, int k) const { return getClampedCell(i,j,k)->getPressure(); }
  
  bool isEmpty(int i, int j, int k) const {
    // for Marker and Cell: does this neighboring cell exist and is it empty?
    if (i < 0 || i >= nx ||
        j < 0 || j >= ny ||
        k < 0 || k >= nz)
      return false;
    return (getCell(i,j,k)->getStatus() == CELL_EMPTY); 
  }
  
  // =================
  // DISPLAY FUNCTIONS
  void Animate();
  void Paint() const; 

  // =================
  // SIMULATION PHASES
  void ComputeNewVelocities();
  float AdjustForIncompressibility();
  void CopyVelocities();
  void MoveParticles();
  void ReassignParticles();
  void SetEmptySurfaceFull();

  // =====================
  // NAVIER-STOKES HELPERS
  float get_u_avg(int i, int j, int k) const { return 0.5*(get_u_plus(i-1,j,k)+get_u_plus(i,j,k)); }
  float get_v_avg(int i, int j, int k) const { return 0.5*(get_v_plus(i,j-1,k)+get_v_plus(i,j,k)); }
  float get_w_avg(int i, int j, int k) const { return 0.5*(get_w_plus(i,j,k-1)+get_w_plus(i,j,k)); }
  float get_u_plus(int i, int j, int k) const { if (u_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_u_plus(); }
  float get_v_plus(int i, int j, int k) const { if (v_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_v_plus(); }
  float get_w_plus(int i, int j, int k) const { if (w_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_w_plus(); }
  float get_new_u_plus(int i, int j, int k) const { if (u_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_new_u_plus(); }
  float get_new_v_plus(int i, int j, int k) const { if (v_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_new_v_plus(); }
  float get_new_w_plus(int i, int j, int k) const { if (w_boundary(i,j,k)) return 0; return getCell(i,j,k)->get_new_w_plus(); }
  void adjust_new_u_plus(int i, int j, int k, float f) { getCell(i,j,k)->adjust_new_u_plus(f); }
  void adjust_new_v_plus(int i, int j, int k, float f) { getCell(i,j,k)->adjust_new_v_plus(f); }
  void adjust_new_w_plus(int i, int j, int k, float f) { getCell(i,j,k)->adjust_new_w_plus(f); }
  float get_uv_plus(int i, int j, int k) const { 
    return 0.5*(get_u_plus(i,j,k) + get_u_plus(i,j+1,k)) * 0.5*(get_v_plus(i,j,k) + get_v_plus(i+1,j,k)); }
  float get_uw_plus(int i, int j, int k) const { 
    return 0.5*(get_u_plus(i,j,k) + get_u_plus(i,j,k+1)) * 0.5*(get_w_plus(i,j,k) + get_w_plus(i+1,j,k)); }
  float get_vw_plus(int i, int j, int k) const { 
    return 0.5*(get_v_plus(i,j,k) + get_v_plus(i,j,k+1)) * 0.5*(get_w_plus(i,j,k) + get_w_plus(i,j+1,k)); }
  bool u_boundary(int i, int j, int k) const { 
    if (i >= 0 && i < nx-1 && j >= 0 && j < ny && k >= 0 && k < nz) return false; return true; }
  bool v_boundary(int i, int j, int k) const { 
    if (i >= 0 && i < nx && j >= 0 && j < ny-1 && k >= 0 && k < nz) return false; return true; }
  bool w_boundary(int i, int j, int k) const { 
    if (i >= 0 && i < nx && j >= 0 && j < ny && k >= 0 && k < nz-1) return false; return true; }

  // ==================
  // DISPLAY HELPER CODE
  void insertColor(float value) const;
  void insertVertex(Vec3f v) const;
  void insertNormal(const Vec3f &p1, const Vec3f &p2, const Vec3f &p3) const;
  void drawTriangle(const Vec3f &p1, const Vec3f &p2, const Vec3f &p3, int wireframe) const;
  void drawIfBoundary(const Vec3f &p1, const Vec3f &p2, const Vec3f &p3, int wireframe) const;

  // ===================
  // MARCHING CUBES CODE
  void MarchingCubes(int x, int y, int z, int wireframe) const;
  void MarchingCubesHelper(Vec3f verts[8], int parity, int wireframe) const;
  void MarchingTetras(Vec3f a, Vec3f b, Vec3f c, Vec3f d, int wireframe) const;
  void MarchingTetrasHelper4(Vec3f a, Vec3f b, Vec3f c, Vec3f d, int wireframe) const;
  void MarchingTetrasHelper3(Vec3f a, Vec3f b, Vec3f c, Vec3f d, int wireframe) const;
  void MarchingTetrasHelper2(Vec3f a, Vec3f b, Vec3f c, Vec3f d, int wireframe) const;
  void MarchingTetrasHelper1(Vec3f a, Vec3f b, Vec3f c, Vec3f d, int wireframe) const;
  float interpolateIsovalue(Vec3f c) const;

  // *********************************************************************  
  // ASSIGNMENT:
  //
  // You can change this function to render isosurfaces of a different scalar field.
  //
  float getIsovalue(int i, int j, int k) const {
    //return getPressure(i,j,k);
    return getClampedCell(i,j,k)->numParticles()/float(density);
  }
  //
  // *********************************************************************  

  float triInterpolate(float x_frac, float y_frac, float z_frac,
		       float a, float b, float c, float d, float e, float f, float g, float h) const;

private:

  // don't use this constructor
  Grid() { assert(0); }
  
  // ==============
  // REPRESENTATION
  ArgParser *args;

  // grid parameters
  int nx,ny,nz;
  float dx;
  float dy;
  float dz;
  Cell *cells;

  // simulation parameters
  Vec3f gravity;
  float dt;
  float viscosity;
  int density; // # of particles we expect to see in a "Full" cell
};

#endif