//#####################################################################
// Copyright 2002-2006, Doug Enright, Ronald Fedkiw, Frederic Gibou, Geoffrey Irving, Frank Losasso, Neil Molino, Avi Robinson-Mosher, Tamar Shinar, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class LEVELSET
//#####################################################################
#ifndef __LEVELSET__
#define __LEVELSET__

#include <Core/Arrays/ARRAYS_FORWARD.h>
#include <Core/Basic_Geometry/BASIC_GEOMETRY_FORWARD.h>
#include <Core/Interpolation/INTERPOLATION_NODE.h>
#include <Core/Interpolation/LINEAR_INTERPOLATION.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/SCALAR_POLICY.h>
#include <cassert>
#include <cfloat>
namespace PhysBAM{

template<class TV> class LEVELSET_CALLBACKS; // TODO: invalid dependency
template<class TV> struct BOUNDARY_POLICY;
template<class TV> struct GRID_ARRAYS_POLICY;
template<class TV> class GRID;
template<class TV,class T2> class BOUNDARY;
template<class TV> class LEVELSET;
template<class T> class INTERVAL;
template<int d> class FACE_INDEX;

template<class TV>
class LEVELSET
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;
public:
    T small_number;
    INTERPOLATION_NODE<TV,T> *interpolation,*curvature_interpolation,*secondary_interpolation;
    INTERPOLATION_NODE<TV,TV>* normal_interpolation;
    T max_time_step;
    T sigma;
    T half_band_width=0;
    bool owns_phi=false;

    BOUNDARY<TV,T>* boundary;
    LEVELSET_CALLBACKS<TV>* levelset_callbacks;
    const ARRAY<bool,FACE_INDEX<TV::m> >* face_velocities_valid_mask_current;
//protected:
    BOUNDARY<TV,T>& boundary_default;
    static LINEAR_INTERPOLATION_NODE<TV,T> interpolation_default;
    static LINEAR_INTERPOLATION_NODE<TV,TV> normal_interpolation_default;
    ARRAY<bool,IV> valid_mask_current;
    ARRAY<bool,IV> valid_mask_next;
private:
    GRID<TV> grid;
public:
    ARRAY<T,IV>& phi;
    ARRAY<TV,IV>* normals;
    ARRAY<T,IV> *curvature;
    ARRAY<INTERVAL<T>,IV> *cell_range;
    int number_of_ghost_cells;

    LEVELSET(int number_of_ghost_cells_input=3);
    LEVELSET(ARRAY<T,IV>& phi_input,const int number_of_ghost_cells_input=3);
    LEVELSET(const GRID<TV>& grid_input,AT_NODES,ARRAY<T,IV>& phi_input,const int number_of_ghost_cells_input=3)
        :LEVELSET(phi_input,number_of_ghost_cells_input)
    {
        Set_Grid(grid_input);
    }
    LEVELSET(const GRID<TV>& grid_input,AT_CELLS,ARRAY<T,IV>& phi_input,const int number_of_ghost_cells_input=3)
        :LEVELSET(phi_input,number_of_ghost_cells_input)
    {
        Set_Grid(grid_input.Get_Regular_Grid_At_MAC_Positions());
    }
    LEVELSET(const LEVELSET&) = delete;
    void operator=(const LEVELSET&) = delete;
    ~LEVELSET();

    const GRID<TV>& Native_Grid() const
    {return grid;}

    void Set_Grid(const GRID<TV>& g)
    {
        grid=g;
        Set_Band_Width();
    }

    const GRID<TV>& Reg_Grid() const
    {
        return grid;
    }
    GRID<TV> Mac_Grid() const
    {
        return grid.Get_MAC_Grid_At_Regular_Positions();
    }
    
    void Set_Small_Number(const T small_number_input=1e-8)
    {small_number=small_number_input;}

    void Set_Max_Time_Step(const T max_time_step_input=1e8)
    {max_time_step=max_time_step_input;}

protected:
    void Initialize_Levelset_Grid_Values(const GRID<TV>& grid_input)
    {}
public:

    void Set_Band_Width(const T number_of_cells=6)
    {half_band_width=number_of_cells*grid.dX.Max()/2;}

    void Initialize_Valid_Masks(const GRID<TV>& grid_input)
    {
        valid_mask_current.Resize(grid_input.Cell_Indices(3),init_new,true);
        valid_mask_next.Resize(grid_input.Cell_Indices(3),no_init);
    }

    void Set_Custom_Boundary(BOUNDARY<TV,T>& boundary_input)
    {boundary=&boundary_input;}

    void Set_Custom_Secondary_Interpolation(INTERPOLATION_NODE<TV,T>& interpolation_input)
    {secondary_interpolation=&interpolation_input;}

    void Set_Custom_Normal_Interpolation(INTERPOLATION_NODE<TV,TV>& normal_interpolation_input)
    {normal_interpolation=&normal_interpolation_input;}

    void Set_Custom_Curvature_Interpolation(INTERPOLATION_NODE<TV,T>& curvature_interpolation_input)
    {curvature_interpolation=&curvature_interpolation_input;}

    void Set_Levelset_Callbacks(LEVELSET_CALLBACKS<TV>& levelset_callbacks_input)
    {levelset_callbacks=&levelset_callbacks_input;}

    void Set_Face_Velocities_Valid_Mask(const ARRAY<bool,FACE_INDEX<TV::m> >* face_velocities_valid_mask_current_input)
    {face_velocities_valid_mask_current=face_velocities_valid_mask_current_input;}

    void Initialize_Levelset_Grid_Values()
    {Initialize_Levelset_Grid_Values(grid);}

    T Phi(const TV& location) const
    {return interpolation->Interpolate_Node(grid,phi,location);}

    T Phi_Secondary(const TV& location) const
    {return secondary_interpolation->Interpolate_Node(grid,phi,location);}

    T Curvature(const TV& location) const // later add finite difference for curvature, like in normal above
    {assert(curvature);return curvature_interpolation->Interpolate_Node(grid,*curvature,location);}

    bool Lazy_Outside(const TV& clamped_X,const T contour_value=0) const
    {return !Lazy_Inside(clamped_X,contour_value);}

    bool Lazy_Outside_Extended_Levelset(const TV& unclamped_X,const T contour_value=0) const
    {return !Lazy_Inside_Extended_Levelset(unclamped_X,contour_value);}

    static TV Normal_At_Node(const GRID<TV>& grid,const ARRAY<T,IV>& phi,const IV& index)
    {TV N;for(int d=0;d<TV::m;d++){IV a(index),b(index);a(d)--;b(d)++;N(d)=(phi(b)-phi(a))*grid.one_over_dX(d);}return N.Normalized();}

    T Compute_Curvature(const IV& index) const
    {return Compute_Curvature(phi,index);}

    RANGE<TV> Bounding_Box() const
    {return grid.domain;}

    T Signed_Distance(const TV& location) const
    {return Extended_Phi(location);}
    
//#####################################################################
    T Extended_Phi(const TV& location) const;
    bool Lazy_Inside(const TV& clamped_X,const T contour_value=0) const;
    bool Lazy_Inside_And_Value(const TV& clamped_X,T& phi_value,const T contour_value=0) const;
    bool Lazy_Inside_Extended_Levelset(const TV& unclamped_X,const T contour_value=0) const;
    // phi_value only guaranteed to be set if function returns true
    bool Lazy_Inside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,const T contour_value=0) const;
    // phi_value only guaranteed to be set if function returns true
    bool Lazy_Outside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,const T contour_value=0) const;
    T CFL(const ARRAY<T,FACE_INDEX<TV::m> >& face_velocities) const;
    T CFL(const ARRAY<TV,IV>& velocity) const;
    TV Iterative_Find_Interface(TV left,TV right,const int iterations=3) const;
    void Compute_Gradient(ARRAY<TV,IV>& gradient,const T time=0) const;
    void Compute_Normals(const T time=0);
    TV Gradient(const TV& location) const;
    TV Normal(const TV& location) const;
    TV Extended_Normal(const TV& location) const;
    SYMMETRIC_MATRIX<T,TV::m> Hessian(const TV& X) const;
    SYMMETRIC_MATRIX<T,TV::m> Hessian(const ARRAY<T,IV>& phi_input,const IV& index) const;
    TV Gradient(const ARRAY<T,IV>& phi_input,const IV& index) const;
    void Compute_Cell_Minimum_And_Maximum(const bool recompute_if_exists=true);
    T Compute_Curvature(const ARRAY<T,IV>& phi_input,const IV& index) const;
    void Compute_Curvature(const T time=0);
    T Compute_Curvature(const TV& location) const;
    void Fast_Marching_Method(const T time=0,const T stopping_distance=0,const ARRAY<IV>* seed_indices=0,const bool add_seed_indices_for_ghost_cells=false,int process_sign=0);
    void Get_Signed_Distance_Using_FMM(ARRAY<T,IV>& signed_distance,const T time=0,const T stopping_distance=0,const ARRAY<IV>* seed_indices=0,
        const bool add_seed_indices_for_ghost_cells=false,int process_sign=0);
    T Approximate_Surface_Size(const T interface_thickness=3,const T time=0) const;
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
//#####################################################################
};
template<class T> class RAY_HIT;
template<class TV> class RAY;

template<class TV> RAY_HIT<typename TV::SCALAR>
Intersects(const RAY<TV>& ray,const LEVELSET<TV>& levelset,typename TV::SCALAR isocontour=0);
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,const LEVELSET<TV>& levelset,const T isocontour=0);
}
#endif
