//#####################################################################
// Copyright 2021.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#ifndef __REFLECTION_HELPER__
#define __REFLECTION_HELPER__

#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAY_VIEW.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/ARRAYS_ND_BASE.h>
#include <Core/Grids/CELL_ITERATOR.h>
#include <Core/Grids/FACE_INDEX.h>
#include <Core/Grids/FACE_ITERATOR.h>
#include <Core/Grids/FACE_RANGE_ITERATOR.h>
#include <Core/Grids/GRID.h>
#include <Core/Grids/NODE_ITERATOR.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Math_Tools/RANGE_ITERATOR.h>
#include <Core/Vectors/VECTOR.h>
#include <functional>

namespace PhysBAM{

template<int d> VECTOR<int,d> Face_To_Fine_Node(const FACE_INDEX<d>& face)
{
    return (2*face.index+1).Add_Axis(face.axis,-1);
}

template<int d> VECTOR<int,d> Cell_To_Fine_Node(const VECTOR<int,d>& cell)
{
    return 2*cell+1;
}

template<int d> VECTOR<int,d> Node_To_Fine_Node(const VECTOR<int,d>& node)
{
    return 2*node;
}

template<int d> FACE_INDEX<d> Fine_Node_To_Face(const VECTOR<int,d>& node,int axis)
{
    return FACE_INDEX<d>(axis,(node.Add_Axis(axis,1)-1)/2);
}

template<int d> VECTOR<int,d> Fine_Node_To_Cell(const VECTOR<int,d>& node)
{
    return (node-1)/2;
}

template<int d> VECTOR<int,d> Fine_Node_To_Node(const VECTOR<int,d>& node)
{
    return node/2;
}

//#####################################################################
// Function Reflection_Cell_Neumann
// Mirrors the cell values in the interior of u_ghost to the ghost region.
//
// If boundary is defined, uses boundary to determine the derivative (normal to the boundary) at a boundary location
// Otherwise assumes that the boundary values are all 0.
//#####################################################################
template<class TV,class T> void Reflection_Cell_Neumann(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Cell_Neumann_Side<TV,T>(grid,u_ghost,ghost,side,boundary);
}
template<class TV,class T> void Reflection_Cell_Neumann_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);
        u_ghost(it.index)=u_ghost(reflected_index);
        if(boundary) u_ghost(it.index)+=boundary(boundary_df,it.side)*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));
    }
}
template<class TV,typename T> void Reflection_Cell_Neumann(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Cell_Neumann_Side<TV,T>(grid,u_ghost,ghost,side,n_boundary);
}
template<class TV,typename T> void Reflection_Cell_Neumann_Side(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost(0).domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost(0).domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);

        for(int a=0;a<TV::m;a++)
        {
            T neumann=n_boundary?n_boundary(boundary_df,it.side,a):(T)0;
            u_ghost(a)(it.index)=u_ghost(a)(reflected_index)+neumann*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));
        }
    }
}

//#####################################################################
// Function Reflection_Cell_NoSlip
// Mirrors and reflects the cell values in the interior of u_ghost to the ghost region around the boundary value.
//
// The boundary values are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//#####################################################################
template<class TV,class T> void Reflection_Cell_NoSlip(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> boundary)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Cell_NoSlip_Side<TV,T>(grid,u_ghost,ghost,side,boundary);
}
template<class TV,class T> void Reflection_Cell_NoSlip_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int)> boundary)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);
        T boundary_value=boundary?boundary(boundary_df,it.side):T();
        u_ghost(it.index)=u_ghost(reflected_index)-2*(u_ghost(reflected_index)-boundary_value);
    }
}
template<class TV,typename T> void Reflection_Cell_NoSlip(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Face_NoSlip_Side<TV,T>(grid,u_ghost,ghost,side,boundary);
}
template<class TV,typename T> void Reflection_Cell_NoSlip_Side(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary)
{
    typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<T,TV::m> TV2;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost(0).domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost(0).domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);

 
        for(int a=0;a<TV2::m;a++)
        {
            T boundary_value=boundary?boundary(boundary_df,it.side,a):(T)0;

            u_ghost(a)(it.index)=u_ghost(a)(reflected_index)-2*(u_ghost(a)(reflected_index)-boundary_value);
        }
    }
}
//#####################################################################
// Function Reflection_Cell_Slip
// Mirrors and reflects the cell values in the interior of u_ghost to the ghost region around the boundary value.
// Assumes that the contents of the cells are vectors of type T.
//    Components of the interior vector that are normal to the boundary are reflected around the boundary value.
//    Components of the interior vector that are tangential to the boundary are mirrored.
//
// The boundary value components are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//    int: the component axis of the boundary value to return
//#####################################################################
template<class TV,class T> void Reflection_Cell_Slip(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int)> n_boundary)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Cell_Slip_Side<TV,T>(grid,u_ghost,ghost,side,d_boundary,n_boundary);
}
template<class TV,class SV> void Reflection_Cell_Slip_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<SV,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<SV(const VECTOR<int,TV::m>&,int)> d_boundary,std::function<SV(const VECTOR<int,TV::m>&,int)> n_boundary)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;       
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);

        // tangential direction
        SV neumann=n_boundary?n_boundary(boundary_df,it.side):SV();
        u_ghost(it.index)=u_ghost(reflected_index)+neumann*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));

        // normal direction
        SV boundary_value=d_boundary?d_boundary(boundary_df,it.side):SV();
        auto value=u_ghost(reflected_index)(axis);
        value-=2*(u_ghost(reflected_index)(axis)-boundary_value(axis));
        u_ghost(it.index)(axis)=value;
    }
}
template<class TV,typename T> void Reflection_Cell_Slip(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Cell_Slip_Side<TV,T>(grid,u_ghost,ghost,side,d_boundary,n_boundary);
}
template<class TV,typename T> void Reflection_Cell_Slip_Side(const GRID<TV>& grid,VECTOR<ARRAY<T,VECTOR<int,TV::m> >,TV::m>& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary)
{
    typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<T,TV::m> TV2;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost(0).domain.Thickened(-ghost));

    for(RANGE_ITERATOR<TV::m> it(u_ghost(0).domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV cell_df=Cell_To_Fine_Node(it.index);
        IV boundary_cell=it.index;
        boundary_cell(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Cell_To_Fine_Node(boundary_cell);
        boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=cell_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Cell(reflected_df);

        for(int a=0;a<TV2::m;a++)
        {
            if(a==axis) // normal direction
            {
                T boundary_value=d_boundary?d_boundary(boundary_df,it.side,a):(T)0;
                T value=u_ghost(a)(reflected_index);
                value-=2*(u_ghost(a)(reflected_index)-boundary_value);
                u_ghost(a)(it.index)=value;
            }
            else // tangential direction
            {
                T neumann=n_boundary?n_boundary(boundary_df,it.side,a):(T)0;
                u_ghost(a)(it.index)=u_ghost(a)(reflected_index)+neumann*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));
            }
        }
    }
}
//#####################################################################
// Function Reflection_Node_Neumann
// Mirrors the node values in the interior of u_ghost to the ghost region.
//
// Assumes that the boundary values are all 0.
//
// Optionally enforces this boundary condition on all interior boundary nodes.
//#####################################################################
template<class TV,class T> void Reflection_Node_Neumann(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Node_Neumann_Side<TV,T>(grid,u_ghost,ghost,side,boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Node_Neumann_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    if(enforce_boundary)
        for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,-ghost,-ghost-1,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
            u_ghost(it.index)=boundary?boundary(Node_To_Fine_Node(it.index),it.side):T();
    
    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV node_df=Node_To_Fine_Node(it.index);
        IV boundary_node=it.index;
        boundary_node(axis)=(it.side&1)?u_size.max_corner(axis)-1:u_size.min_corner(axis);
        IV boundary_df=Node_To_Fine_Node(boundary_node);
        IV reflected_df=boundary_df;
        reflected_df(axis)-=node_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Node(reflected_df);
        u_ghost(it.index)=u_ghost(reflected_index);
        if(boundary) u_ghost(it.index)+=boundary(boundary_df,it.side)*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));
    }
}
//#####################################################################
// Function Reflection_Node_NoSlip
// Mirrors and reflects the node values in the interior of u_ghost to the ghost region around the boundary value.
//
// The boundary values are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//
// Optionally enforces the boundary values provided on the interior boundary nodes.
//#####################################################################
template<class TV,class T> void Reflection_Node_NoSlip(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Node_NoSlip_Side<TV,T>(grid,u_ghost,ghost,side,boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Node_NoSlip_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int)> boundary=0,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    if(enforce_boundary)
        for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,-ghost,-ghost-1,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
            u_ghost(it.index)=boundary?boundary(Node_To_Fine_Node(it.index),it.side):T();

    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV node_df=Node_To_Fine_Node(it.index);
        IV boundary_node=it.index;
        boundary_node(axis)=(it.side&1)?u_size.max_corner(axis)-1:u_size.min_corner(axis);
        IV boundary_df=Node_To_Fine_Node(boundary_node);
        IV reflected_df=boundary_df;
        reflected_df(axis)-=node_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Node(reflected_df);
        T boundary_value=boundary?boundary(boundary_df,it.side):T();
        u_ghost(it.index)=u_ghost(reflected_index)-2*(u_ghost(reflected_index)-boundary_value);
    }
}
//#####################################################################
// Function Reflection_Node_Slip
// Mirrors and reflects the node values in the interior of u_ghost to the ghost region around the boundary value.
// Assumes that the contents of the nodes are vectors of type T.
//    Components of the interior vector that are normal to the boundary are reflected around the boundary value.
//    Components of the interior vector that are tangential to the boundary are mirrored.
//
// The boundary value components are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//    int: the component axis of the boundary value to return
//
// Optionally enforces the boundary values provided on the interior boundary nodes.
//#####################################################################
template<class TV,class T> void Reflection_Node_Slip(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int)> n_boundary,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Node_Slip_Side<TV,T>(grid,u_ghost,ghost,side,d_boundary,n_boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Node_Slip_Side(const GRID<TV>& grid,ARRAYS_ND_BASE<T,VECTOR<int,TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int)> n_boundary,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<T,TV::m> TV2;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    if(enforce_boundary)
        for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,-ghost,-ghost-1,RI::ghost,side);it.Valid();it.Next())
            u_ghost(it.index)=d_boundary?d_boundary(Node_To_Fine_Node(it.index),it.side):T();

    for(RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RI::ghost|RI::delay_corners,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        IV node_df=Node_To_Fine_Node(it.index);
        IV boundary_node=it.index;
        boundary_node(axis)=(it.side&1)?u_size.max_corner(axis)-1:u_size.min_corner(axis);
        IV boundary_df=Node_To_Fine_Node(boundary_node);
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=node_df(axis)-boundary_df(axis);
        IV reflected_index=Fine_Node_To_Node(reflected_df);

        // tangential direction
        T neumann=n_boundary?n_boundary(boundary_df,it.side):T();
        u_ghost(it.index)=u_ghost(reflected_index)+neumann*grid.dX(axis)*abs(reflected_index(axis)-it.index(axis));

        // normal direction
        T boundary_value=d_boundary?d_boundary(boundary_df,it.side):T();
        auto value=u_ghost(reflected_index)(axis);
        value-=2*(u_ghost(reflected_index)(axis)-boundary_value(axis));
        u_ghost(it.index)(axis)=value;
    }
}
//#####################################################################
// Function Reflection_Face_Neumann
// Mirrors the face values in the interior of u_ghost to the ghost region.
//
// If boundary is defined, uses boundary to determine the derivative (normal to the boundary) at a boundary location
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//    int: the axis of the face index 
// Otherwise assumes that the boundary derivative values are all 0.
//#####################################################################
template<class TV,class T> void Reflection_Face_Neumann(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary=0)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Face_Neumann_Side<TV,T>(grid,u_ghost,ghost,side,boundary);
}
template<class TV,class T> void Reflection_Face_Neumann_Side(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary=0)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));
    
    for(FACE_RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RF::ghost|RF::delay_corners|RF::skip_inner,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        int face_axis=it.face.axis;
        IV face_df=Face_To_Fine_Node(it.face);
        FACE_INDEX<TV::m> boundary_face=it.face;
        boundary_face.index(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Face_To_Fine_Node(boundary_face);
        if(face_axis!=axis) boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=face_df(axis)-boundary_df(axis);
        FACE_INDEX<TV::m> reflected_index=Fine_Node_To_Face(reflected_df,face_axis);
        u_ghost(it.face)=u_ghost(reflected_index);
        if(boundary) u_ghost(it.face)+=boundary(boundary_df,it.side,face_axis)*grid.dX(axis)*abs(reflected_index.index(axis)-it.face.index(axis));
    }
}
//#####################################################################
// Function Reflection_Face_NoSlip
// Mirrors and reflects the face values in the interior of u_ghost to the ghost region around the boundary value.
//
// The boundary values are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//    int: the axis of the face index 
//
// Optionally enforces the boundary values provided on the interior boundary faces.
//#####################################################################
template<class TV,class T> void Reflection_Face_NoSlip(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Face_NoSlip_Side<TV,T>(grid,u_ghost,ghost,side,boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Face_NoSlip_Side(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    if(enforce_boundary)
        for(FACE_RANGE_ITERATOR<TV::m> it(u_size,0,0,RF::ghost,side);it.Valid();it.Next())
            u_ghost(it.face)=boundary?boundary(Face_To_Fine_Node(it.face),it.side,it.face.axis):T();

    for(FACE_RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RF::ghost|RF::delay_corners|RF::skip_inner,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        int face_axis=it.face.axis;
        IV face_df=Face_To_Fine_Node(it.face);
        FACE_INDEX<TV::m> boundary_face=it.face;
        boundary_face.index(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Face_To_Fine_Node(boundary_face);
        if(face_axis!=axis) boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=face_df(axis)-boundary_df(axis);
        FACE_INDEX<TV::m> reflected_index=Fine_Node_To_Face(reflected_df,face_axis);
        T boundary_value=boundary?boundary(boundary_df,it.side,face_axis):T();

        u_ghost(it.face)=u_ghost(reflected_index)-2*(u_ghost(reflected_index)-boundary_value);
    }
}
//#####################################################################
// Function Reflection_Face_Slip
// Mirrors and reflects the face values in the interior of u_ghost to the ghost region around the boundary value.
//    Values of the interior faces that are normal to the boundary are reflected around the boundary value.
//    Values of the interior faces that are tangential to the boundary are mirrored.
//
// The boundary values are provided by the user-provided function 'boundary'
// boundary parameters:
//    VECTOR<int,m>: boundary location index (into a double-fine version of the provided grid)
//    int: the side of the grid that the boundary location is on
//    int: the axis of the face index 
//
// Optionally enforces the boundary values provided on the interior boundary faces.
//#####################################################################
template<class TV,class T> void Reflection_Face_Slip(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Face_Slip_Side<TV,T>(grid,u_ghost,ghost,side,boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Face_Slip_Side(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> boundary,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));

    if(enforce_boundary)
        for(FACE_RANGE_ITERATOR<TV::m> it(u_size,0,0,RF::ghost,side);it.Valid();it.Next())
            u_ghost(it.face)=boundary?boundary(Face_To_Fine_Node(it.face),it.side,it.face.axis):T();

    for(FACE_RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RF::ghost|RF::delay_corners|RF::skip_inner,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        int face_axis=it.face.axis;
        IV face_df=Face_To_Fine_Node(it.face);
        FACE_INDEX<TV::m> boundary_face=it.face;
        boundary_face.index(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Face_To_Fine_Node(boundary_face);
        if(face_axis!=axis) boundary_df(axis)-=1;

        IV reflected_df=boundary_df;
        reflected_df(axis)-=face_df(axis)-boundary_df(axis);
        FACE_INDEX<TV::m> reflected_index=Fine_Node_To_Face(reflected_df,face_axis);
        T boundary_value=boundary?boundary(boundary_df,it.side,face_axis):T();

        T value=u_ghost(reflected_index);
        if(axis==face_axis) value-=2*(u_ghost(reflected_index)-boundary_value);
        u_ghost(it.face)=value;
    }
}
template<class TV,class T> void Reflection_Face_Slip(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,std::function<T(const VECTOR<int,TV::m>&,int,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary,bool enforce_boundary=false)
{
    for(int side=0;side<TV::m*2;side++)
        Reflection_Face_Slip_Side<TV,T>(grid,u_ghost,ghost,side,d_boundary,n_boundary,enforce_boundary);
}
template<class TV,class T> void Reflection_Face_Slip_Side(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& u_ghost,int ghost,int side,std::function<T(const VECTOR<int,TV::m>&,int,int)> d_boundary,std::function<T(const VECTOR<int,TV::m>&,int,int)> n_boundary,bool enforce_boundary=false)
{
    typedef VECTOR<int,TV::m> IV;
    GRID<TV> grid_df(grid.numbers_of_cells*2,grid.domain);
    RANGE<IV> u_size(u_ghost.domain.Thickened(-ghost));
    if(enforce_boundary)
        for(FACE_RANGE_ITERATOR<TV::m> it(u_size,0,0,RF::ghost,side);it.Valid();it.Next())
            u_ghost(it.face)=d_boundary?d_boundary(Face_To_Fine_Node(it.face),it.side,it.face.axis):T();

    for(FACE_RANGE_ITERATOR<TV::m> it(u_ghost.domain,0,-ghost,RF::ghost|RF::delay_corners|RF::skip_inner,side);it.Valid();it.Next())
    {
        int axis=it.side/2;
        int face_axis=it.face.axis;
        IV face_df=Face_To_Fine_Node(it.face);
        FACE_INDEX<TV::m> boundary_face=it.face;
        boundary_face.index(axis)=(it.side&1)?u_size.max_corner(axis):u_size.min_corner(axis);
        IV boundary_df=Face_To_Fine_Node(boundary_face);
        if(face_axis!=axis) boundary_df(axis)-=1;
        
        IV reflected_df=boundary_df;
        reflected_df(axis)-=face_df(axis)-boundary_df(axis);
        FACE_INDEX<TV::m> reflected_index=Fine_Node_To_Face(reflected_df,face_axis);
        if(axis==face_axis) // normal direction
        {
            T boundary_value=d_boundary?d_boundary(boundary_df,it.side,face_axis):T();
            T value=u_ghost(reflected_index);
            value-=2*(u_ghost(reflected_index)-boundary_value);
            u_ghost(it.face)=value;
        }
        else // tangential direction
        {
            T neumann=n_boundary?n_boundary(boundary_df,it.side,face_axis):T();
            u_ghost(it.face)=u_ghost(reflected_index)+neumann*grid.dX(axis)*abs(reflected_index.index(axis)-it.face.index(axis));
        }
    }
}
}
//#####################################################################
#endif
