//#####################################################################
// Copyright 2023, Mehrnaz Ayazi.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################

#ifndef __WENO_INTERPOLATION_FACE__
#define __WENO_INTERPOLATION_FACE__

#include <Core/Grids/FACE_INDEX.h>
#include <Core/Interpolation/DIV_FREE_BALSARA.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>
#include <Core/Interpolation/INTERPOLATION_HELPERS.h>
#include <Core/Interpolation/WENO_INTERPOLATION.h>
#include <Core/Vectors/VECTOR.h>

namespace PhysBAM
{
template<class TV,class T2,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
class WENO_INTERPOLATION_FACE:public INTERPOLATION_FACE<TV,T2>
{
    typedef typename TV::SCALAR T;typedef VECTOR<int, TV::m> IV;
    typedef VECTOR<T2,TV::m> TV2;typedef SYMMETRIC_MATRIX<T2,TV::m> SM2;
    typedef MATRIX<T2, TV::m> TM2;typedef DIAGONAL_MATRIX<T2,TV::m> DM2;
    typedef INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY> BASE;
    using typename INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY>::REF;

public:
    T eps=(T)1e-6;
    
    WENO_INTERPOLATION_FACE() = default;

    TV2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {
        REF ref(u,X);
        auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<2,false>(grid,u.Domain_Indices(),X);
        TV2 r;
        for(int i=0;i<TV::m;i++)
            r(i)=WENO_Interpolation<T,T2,TV::m>(p(i).second,u.Component(i),p(i).first,eps);
        return r;
    }

    T2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {
        REF ref(u,X);
        auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<2,false>(grid,u.Domain_Indices(),X,axis);
        return WENO_Interpolation<T,T2,TV::m>(p.second,u.Component(axis),p.first,eps);
    }

    VECTOR<VECTOR<T2,2>,TV::m> Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X) const
    {PHYSBAM_NOT_IMPLEMENTED();}
    VECTOR<T2,2> Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X,int axis) const
    {PHYSBAM_NOT_IMPLEMENTED();}

    VECTOR<TV2,TV::m> Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {PHYSBAM_NOT_IMPLEMENTED();}
    TV2 Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {PHYSBAM_NOT_IMPLEMENTED();}

    VECTOR<SM2,TV::m> Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {PHYSBAM_NOT_IMPLEMENTED();}
    SM2 Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {PHYSBAM_NOT_IMPLEMENTED();}

protected:
    virtual TV2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const override
    {return Interpolate_Face(grid,u,X);}
    virtual T2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const override
    {return Interpolate_Face(grid,u,X,axis);}
    virtual VECTOR<VECTOR<T2,2>,TV::m> Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,
        const T_FACE_ARRAY& u1,const TV& X) const override
    {return Extrema_Face(grid,u0,u1,X);}
    virtual VECTOR<T2,2> Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,
        const T_FACE_ARRAY& u1,const TV& X,int axis) const override
    {return Extrema_Face(grid,u0,u1,X,axis);}
    virtual VECTOR<TV2,TV::m> Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const override
    {return Gradient_Face(grid,u,X);}
    virtual TV2 Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const override
    {return Gradient_Face(grid,u,X,axis);}
    virtual VECTOR<SM2,TV::m> Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const override
    {return Hessian_Face(grid,u,X);}
    virtual SM2 Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const override
    {return Hessian_Face(grid,u,X,axis);}
//#####################################################################
};
}
#endif
