//#####################################################################
// Copyright 2002-2008, Zhaosheng Bao, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Sergey Koltakov, Nipun Kwatra, Frank Losasso, Andrew Selle, Tamar Shinar, Eftychios Sifakis, Jerry Talton, Joseph Teran.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class GRID
//#####################################################################
// A number of functions (e.g. the Clamp functions) assume the grid indexing starts at 1.  This way we can use truncation rather than floor because floor is really slow.
//#####################################################################
#ifndef __GRID__
#define __GRID__

#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Grids/FACE_INDEX.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{

template<class T> class GRID_ITERATOR_BASE;
template<class T> class NODE_ITERATOR;
template<class T> class CELL_ITERATOR;
template<class T> class FACE_ITERATOR;

template<bool is_mac> struct GRID_TYPE{};
typedef GRID_TYPE<0> AT_NODES;
typedef GRID_TYPE<1> AT_CELLS;
static const GRID_TYPE<0> reg_grid_type,at_nodes;
static const GRID_TYPE<1> mac_grid_type,at_cells;

enum GRID_REGION {WHOLE_REGION,GHOST_REGION,BOUNDARY_REGION,INTERIOR_REGION,BOUNDARY_INTERIOR_REGION}; // for iterators

template<class TV>
class GRID
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;
    friend class GRID<VECTOR<T,1> >;
    friend class GRID<VECTOR<T,2> >;
    friend class GRID<VECTOR<T,3> >;
public:
    static const int dimension=TV::m;
    static const int number_of_cells_per_block=1<<dimension;
    static const int number_of_faces_per_block=dimension==1?3:dimension==2?12:36;
    static const int number_of_incident_faces_per_block=2*dimension-2;
    static const int number_of_nodes_per_face=1<<(dimension-(dimension!=0));
    static const int number_of_nodes_per_cell=1<<dimension;
    static const int number_of_cells_per_node=1<<dimension;
    static const int number_of_neighbors_per_node=2*dimension;
    static const int number_of_neighbors_per_cell=2*dimension;
    static const int number_of_faces_per_cell=2*dimension;
    static const int number_of_one_ring_neighbors_per_cell=dimension==1?2:dimension==2?8:26;

    typedef T SCALAR;
    typedef TV VECTOR_T;
    typedef IV VECTOR_INT;
    typedef IV INDEX;

    RANGE<TV> domain;
    TV dX;
    TV one_over_dX;
    IV numbers_of_cells; // not saved to file

    GRID()
        :domain(RANGE<TV>::Unit_Box())
    {
    }

    GRID(const IV& cells,const RANGE<TV>& box)
    {
        Initialize(cells,box);
    }

private:
    GRID(const IV& cells,const TV& dX,const TV& one_over_dX,const RANGE<TV>& box)
        :domain(box),dX(dX),one_over_dX(one_over_dX),numbers_of_cells(cells)
    {
    }
public:

    template<class T2>
    GRID(const GRID<T2>& grid_input)
    {
        Initialize(grid_input.numbers_of_cells,RANGE<TV>(grid_input.domain));
    }

    GRID<TV>& operator=(const GRID<TV>& grid_input) = default;

    bool operator==(const GRID<TV>& grid) const
    {return grid.numbers_of_cells==numbers_of_cells && grid.domain==domain;}

    inline bool operator!=(const GRID<TV>& grid) const
    {return !operator==(grid);}

    const TV& DX() const
    {return dX;}

    const TV& One_Over_DX() const
    {return one_over_dX;}

    const IV& Numbers_Of_Cells() const
    {return numbers_of_cells;}

    IV Numbers_Of_Nodes() const
    {return numbers_of_cells+1;}

    T Cell_Size() const
    {return dX.Product();}

    TV Face_Sizes() const
    {return Cell_Size()*one_over_dX;}

    T Face_Size(const int axis) const
    {return Cell_Size()*one_over_dX(axis);}

    // Location of center of cell
    TV Cell_X(const IV& index) const
    {return domain.min_corner+(TV(index)+(T).5)*dX;}

    // Index of nearest cell center
    IV Cell_I(const TV& X) const
    {return Node_L(X);}

    // Index of cell center to the lower left of X
    IV Cell_L(const TV& X) const
    {return IV(floor((X-domain.min_corner)*one_over_dX-(T).5));}

    // Index of cell center to the upper right of X
    IV Cell_U(const TV& X) const
    {return IV(ceil((X-domain.min_corner)*one_over_dX-(T).5));}

    // Index of nearest cell center, clamped to valid indices such that
    // index+inset_max is also a valid index in a grid with ghost
    // cells. (inset_max>=0)
    IV Cell_IC(const TV& X,int ghost=0,bool inset_max=false) const
    {return clamp(Cell_I(X),IV()-ghost,numbers_of_cells+(ghost-inset_max-1));}

    // Index of lower left cell center, clamped to valid indices such that
    // index+inset_max is also a valid index in a grid with ghost
    // cells. (inset_max>=0)
    IV Cell_LC(const TV& X,int ghost=0,bool inset_max=false) const
    {return clamp(Cell_L(X),IV()-ghost,numbers_of_cells+(ghost-inset_max-1));}

    // Location of node
    TV Node_X(const IV& index) const
    {return domain.min_corner+TV(index)*dX;}

    // Index of nearest node
    IV Node_I(const TV& X) const
    {return IV(floor((X-domain.min_corner)*one_over_dX+(T).5));}

    // Index of node center to the upper left of X
    IV Node_L(const TV& X) const
    {return IV(floor((X-domain.min_corner)*one_over_dX));}

    // Index of node center to the lower right of X
    IV Node_U(const TV& X) const
    {return IV(ceil((X-domain.min_corner)*one_over_dX));}

    // Index of nearest node, clamped to valid indices such that
    // index+inset_max is also a valid index in a grid with ghost
    // cells. (inset_max>=0)
    IV Node_IC(const TV& X,int ghost=0,bool inset_max=false) const
    {return clamp(Node_I(X),IV()-ghost,numbers_of_cells+(ghost-inset_max));}

    // Index of lower left node, clamped to valid indices such that
    // index+inset_max is also a valid index in a grid with ghost
    // cells. (inset_max>=0)
    IV Node_LC(const TV& X,int ghost=0,bool inset_max=false) const
    {return clamp(Node_L(X),IV()-ghost,numbers_of_cells+(ghost-inset_max));}

    // Location of face center
    TV Face_X(const FACE_INDEX<TV::m>& face) const
    {return Cell_X(face.index).Add_Axis(face.axis,-(T).5*dX(face.axis));}

    // Index of nearest face
    FACE_INDEX<TV::m> Face_I(const TV& X,int axis) const
    {return FACE_INDEX<TV::m>(axis,Cell_I(X.Add_Axis(axis,(T).5*dX(axis))));}

    // range of cells that entirely contains input range; half-open.
    RANGE<IV> Cell_Range(const RANGE<TV>& range) const
    {return RANGE<IV>(Cell_I(range.min_corner),Cell_I(range.max_corner)+1);}

    RANGE<TV> Face_Domain(const int axis,const IV& index,T thickness_over_two=0) const
    {TV dimensions=(T).5*dX;dimensions(axis)=thickness_over_two;RANGE<TV> domain(Face(axis,index));domain.Change_Size(dimensions);return domain;}

    static IV Node_Cell_Index(const IV& node,int cell)
    {return node+Binary_Counts(node)(cell)-1;}

    static IV Face_Node_Index(const int axis,const VECTOR<int,3>& face_index,int node)
    {static const IV corner_from_face_offset[3][4]={
        {IV(0,0,0),IV(0,1,0),IV(0,0,1),IV(0,1,1)},
        {IV(0,0,0),IV(1,0,0),IV(0,0,1),IV(1,0,1)},
        {IV(0,0,0),IV(1,0,0),IV(0,1,0),IV(1,1,0)}};
    assert((unsigned)axis<3&&(unsigned)node<4);return face_index+corner_from_face_offset[axis][node];}

    static IV Face_Node_Index(const int axis,const VECTOR<int,2>& face_index,int node)
    {assert((unsigned)node<2);IV index=face_index;index[1-axis]+=node;return index;}

    static IV Face_Node_Index(const int axis,const VECTOR<int,1>& face_index,int node)
    {assert(axis==0&&node==0);return face_index;}

    static IV Node_Face_Index(const int axis,const VECTOR<int,3>& node_index,int face)
    {static const IV face_from_node_offset[3][4]={
        {IV(0,-1,-1),IV(0,0,-1),IV(0,-1,0),IV(0,0,0)},
        {IV(-1,0,-1),IV(-1,0,0),IV(0,0,-1),IV(0,0,0)},
        {IV(-1,-1,0),IV(0,-1,0),IV(-1,0,0),IV(0,0,0)}};
    assert((unsigned)face<4&&(unsigned)axis<3);return node_index+face_from_node_offset[axis][face];}

    static IV Node_Face_Index(const int axis,const VECTOR<int,2>& node_index,int face)
    {assert((unsigned)face<2);IV index=node_index;index[1-axis]+=face-1;return index;}

    static IV Node_Face_Index(const int axis,const VECTOR<int,1>& node_index,int face)
    {assert(axis==0&&face==0);return node_index;}

    RANGE<TV> Cell_Domain(const IV& index) const
    {TV corner=Node_X(index);return RANGE<TV>(corner,corner+dX);}

    RANGE<TV> Cell_Domain(const RANGE<IV>& cell_range) const
    {return RANGE<TV>(Node_X(cell_range.min_corner),Node_X(cell_range.max_corner));}

    static void Cells_Touching_Face(const int axis,const IV& face_index,IV& cell1,IV& cell2)
    {cell2=face_index;cell1=face_index;cell1(axis)--;}

    RANGE<IV> Clamp_To_Cell(const RANGE<TV>& box,int ghost=0) const
    {return RANGE<IV>(Cell_IC(box.min_corner,ghost),Cell_IC(box.max_corner,ghost)+1);}

    const RANGE<TV>& Domain() const
    {return domain;}

    RANGE<TV> Ghost_Domain(const int ghost) const
    {return domain.Thickened(dX*ghost);}

    RANGE<IV> Node_Indices(int ghost=0) const
    {return RANGE<IV>(IV(),Numbers_Of_Nodes()).Thickened(ghost);}

    RANGE<IV> Cell_Indices(int ghost=0) const
    {return RANGE<IV>(IV(),numbers_of_cells).Thickened(ghost);}

    RANGE<IV> Face_Indices(int axis,int ghost=0) const
    {auto ci=Cell_Indices(ghost);ci.max_corner(axis)++;return ci;}

    bool Valid_Cell(const IV& index,int ghost=0) const
    {return Cell_Indices(ghost).Lazy_Inside_Half_Open(index);}

    bool Valid_Node(const IV& index,int ghost=0) const
    {return Node_Indices(ghost).Lazy_Inside_Half_Open(index);}

    bool Valid_Face(const FACE_INDEX<TV::m>& face,int ghost=0) const
    {return Face_Indices(face.axis,ghost).Lazy_Inside_Half_Open(face.index);}

    bool Outside(const TV& X) const
    {return domain.Lazy_Outside(X);}

    GRID<TV> Get_Face_Grid(const int axis) const
    {IV numbers=numbers_of_cells;numbers(axis)++;TV offset((T).5*dX);offset(axis)=0;
    return {IV().Componentwise_Max(numbers)-1,dX,one_over_dX,RANGE<TV>(domain.min_corner+offset,domain.max_corner-offset)};}

    GRID<TV> Get_Face_MAC_Grid(const int axis) const
    {IV numbers=numbers_of_cells;numbers(axis)++;TV offset;offset(axis)=(T).5*dX(axis);
    return {IV().Componentwise_Max(numbers),dX,one_over_dX,RANGE<TV>(domain.min_corner-offset,domain.max_corner+offset)};}
    
    GRID<TV> Get_Center_Grid() const
    {return {numbers_of_cells-1,dX,one_over_dX,domain.Thickened(-(T).5*dX)};}

    GRID<TV> Get_Regular_Grid_At_MAC_Positions() const
    {return Get_Center_Grid();}

    GRID<TV> Get_MAC_Grid_At_Regular_Positions() const
    {return {numbers_of_cells+1,dX,one_over_dX,domain.Thickened((T).5*dX)};}

    GRID<TV> Get_Ghost_Grid(int ghost) const
    {return {numbers_of_cells+2*ghost,dX,one_over_dX,domain.Thickened(dX*ghost)};}
    
    GRID<VECTOR<T,TV::m-1>> Remove_Dimension(int dimension) const
    {return {numbers_of_cells.Remove_Index(dimension),dX.Remove_Index(dimension),
         one_over_dX.Remove_Index(dimension),domain.Remove_Dimension(dimension)};}

    GRID<VECTOR<T,1>> Get_1D_Grid(int axis) const
    {return {VECTOR<int,1>(numbers_of_cells(axis)),VECTOR<T,1>(dX(axis)),VECTOR<T,1>(one_over_dX(axis)),
         RANGE<VECTOR<T,1> >(VECTOR<T,1>(domain.min_corner(axis)),VECTOR<T,1>(domain.max_corner(axis)))};}

    static VECTOR<int,3> Node_Neighbor(const VECTOR<int,3>& index,int i)
    {static const VECTOR<int,3> neighbor_offset[6]={VECTOR<int,3>(-1,0,0),VECTOR<int,3>(1,0,0),VECTOR<int,3>(0,-1,0),VECTOR<int,3>(0,1,0),VECTOR<int,3>(0,0,-1),VECTOR<int,3>(0,0,1)};
    assert((unsigned)i<6);return index+neighbor_offset[i];}

    static VECTOR<int,2> Node_Neighbor(const VECTOR<int,2>& index,int i)
    {static const VECTOR<int,2> neighbor_offset[4]={VECTOR<int,2>(-1,0),VECTOR<int,2>(1,0),VECTOR<int,2>(0,-1),VECTOR<int,2>(0,1)};
    assert((unsigned)i<4);return index+neighbor_offset[i];}

    static VECTOR<int,1> Node_Neighbor(const VECTOR<int,1>& index,int i) // i=0 to 1
    {static const VECTOR<int,1> neighbor_offset[2]={VECTOR<int,1>(-1),VECTOR<int,1>(1)};
    assert((unsigned)i<2);return index+neighbor_offset[i];}

    static VECTOR<int,0> Node_Neighbor(const VECTOR<int,0>& index,int i) // i=0
    {assert(0==i);return index;}

    static void Cells_Neighboring_Node(const VECTOR<int,3>& node,VECTOR<int,3> cells[8])
    {cells[0]=VECTOR<int,3>(node.x-1,node.y-1,node.z-1);cells[1]=VECTOR<int,3>(node.x,node.y-1,node.z-1);cells[2]=VECTOR<int,3>(node.x-1,node.y,node.z-1);
    cells[3]=VECTOR<int,3>(node.x,node.y,node.z-1);cells[4]=VECTOR<int,3>(node.x-1,node.y-1,node.z);cells[5]=VECTOR<int,3>(node.x,node.y-1,node.z);
    cells[6]=VECTOR<int,3>(node.x-1,node.y,node.z);cells[7]=VECTOR<int,3>(node.x,node.y,node.z);}

    static void Cells_Neighboring_Node(const VECTOR<int,2>& node,VECTOR<int,2> cells[4])
    {cells[0]=IV(node.x-1,node.y-1);cells[1]=IV(node.x,node.y-1);cells[2]=IV(node.x-1,node.y);cells[3]=IV(node.x,node.y);}

    static void Cells_Neighboring_Node(const VECTOR<int,1>& node,VECTOR<int,1> cells[2])
    {cells[0]=IV(node.x-1);cells[1]=IV(node.x);}

    static void Nodes_In_Cell_From_Minimum_Corner_Node(const VECTOR<int,3>& minimum_corner_node,VECTOR<int,3> nodes[8])
    {int i=minimum_corner_node.x,j=minimum_corner_node.y,k=minimum_corner_node.z;
    nodes[0]=minimum_corner_node;nodes[1]=VECTOR<int,3>(i+1,j,k);nodes[2]=VECTOR<int,3>(i,j+1,k);nodes[3]=VECTOR<int,3>(i+1,j+1,k);
    nodes[4]=VECTOR<int,3>(i,j,k+1);nodes[5]=VECTOR<int,3>(i+1,j,k+1);nodes[6]=VECTOR<int,3>(i,j+1,k+1);nodes[7]=VECTOR<int,3>(i+1,j+1,k+1);}

    static void Nodes_In_Cell_From_Minimum_Corner_Node(const VECTOR<int,2>& minimum_corner_node,VECTOR<int,2> nodes[4])
    {nodes[0]=minimum_corner_node;nodes[1]=VECTOR<int,2>(minimum_corner_node.x+1,minimum_corner_node.y);nodes[2]=VECTOR<int,2>(minimum_corner_node.x,minimum_corner_node.y+1);
    nodes[3]=VECTOR<int,2>(minimum_corner_node.x+1,minimum_corner_node.y+1);}

    static void Nodes_In_Cell_From_Minimum_Corner_Node(const VECTOR<int,1>& minimum_corner_node,VECTOR<int,1> nodes[2])
    {nodes[0]=minimum_corner_node;nodes[1]=VECTOR<int,1>(minimum_corner_node.x+1);}

    static IV First_Face_Index_In_Cell(const int axis,const IV& cell_index)
    {return cell_index;}

    static IV Second_Face_Index_In_Cell(const int axis,const IV& cell_index)
    {IV face_index=cell_index;face_index(axis)++;return face_index;}

    static IV One_Ring_Neighbor(const VECTOR<int,3>& index,int i)
    {static const IV neighbor_offset[26]={IV(-1,-1,-1),IV(0,-1,-1),IV(1,-1,-1),IV(-1,0,-1),IV(0,0,-1),
    IV(1,0,-1),IV(-1,1,-1),IV(0,1,-1),IV(1,1,-1),IV(-1,-1,0),IV(0,-1,0),IV(1,-1,0),
    IV(-1,0,0),IV(1,0,0),IV(-1,1,0),IV(0,1,0),IV(1,1,0),IV(-1,-1,1),IV(0,-1,1),
    IV(1,-1,1),IV(-1,0,1),IV(0,0,1),IV(1,0,1),IV(-1,1,1),IV(0,1,1), IV(1,1,1)};
    assert((unsigned)i<26);return index+neighbor_offset[i];}

    static IV One_Ring_Neighbor(const VECTOR<int,2>& index,int i)
    {static const IV neighbor_offset[8]={IV(-1,-1),IV(0,-1),IV(1,-1),IV(-1,0),IV(1,0),IV(-1,1),
        IV(0,1),IV(1,1)};
    assert((unsigned)i<8);return index+neighbor_offset[i];}

    static IV One_Ring_Neighbor(const VECTOR<int,1>& index,int i)
    {return Node_Neighbor(index,i);}

    static void Neighboring_Faces(VECTOR<FACE_INDEX<TV::m>,number_of_faces_per_cell>& n,const IV& index)
    {for(int a=0;a<TV::m;a++){FACE_INDEX<TV::m> fi(a,index);n(a*2)=fi;fi.index(a)++;n(a*2+1)=fi;}}

//#####################################################################
    void Initialize(const IV& cells_input,const RANGE<TV>& box);
    static GRID<TV> Create_Grid_Given_Cell_Size(const RANGE<TV>& domain,T cell_size,int boundary_nodes=0);
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
//#####################################################################
};
template<class TV>
inline std::ostream& operator<<(std::ostream& output,const GRID<TV>& grid)
{output<<"("<<grid.domain<<" "<<grid.numbers_of_cells<<")";return output;}
}
#endif
