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

#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>
namespace PhysBAM{

template<class TV,class T2,class STENCIL,bool const_extrap=false,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
class GENERAL_INTERPOLATION_FACE:public INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY>
{
    static_assert(TV::m>=1 && TV::m<=3);
    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 INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY> BASE;
    constexpr static int ss=STENCIL::stencil_size;
    using typename INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY>::REF;

  public:
    typedef T_FACE_ARRAY FACE_LOOKUP;
    STENCIL stencil;
    
    GENERAL_INTERPOLATION_FACE() = default;
    GENERAL_INTERPOLATION_FACE(const GENERAL_INTERPOLATION_FACE&) = delete;
    void operator=(const GENERAL_INTERPOLATION_FACE&) = delete;
    virtual ~GENERAL_INTERPOLATION_FACE() = default;

    TV2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    T2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;

    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;
    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;

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

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

protected:
    virtual TV2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual T2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
    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;
    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;
    virtual VECTOR<TV2,TV::m> Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual TV2 Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
    virtual VECTOR<SM2,TV::m> Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const;
    virtual SM2 Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const;
//#####################################################################
};
}
#endif
