//#####################################################################
// 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 __INTERPOLATION_FACE__
#define __INTERPOLATION_FACE__

#include <Core/Arrays/ARRAY.h>
#include <Core/Grids/GRID.h>
#include <Core/Interpolation/FACE_LOOKUP_UNIFORM.h>
#include <Core/Log/DEBUG_UTILITIES.h>
namespace PhysBAM{

template<class TV,class T2,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
class INTERPOLATION_FACE
{
    typedef typename TV::SCALAR T;typedef VECTOR<int,TV::m> IV;
    typedef VECTOR<T2,TV::m> TV2;typedef SYMMETRIC_MATRIX<T2,TV::m> SM2;

public:
    INTERPOLATION_FACE() = default;
    INTERPOLATION_FACE(const INTERPOLATION_FACE&) = delete;
    void operator=(const INTERPOLATION_FACE&) = delete;
    virtual ~INTERPOLATION_FACE() = default;
    
    // These are overridden by each base class to allow compile-time function
    // calls to implementations.  The base class implementation calls the
    // virtual versions below to do dynamic lookup when the base type is not
    // known.
    TV2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {return Interpolate_Face_V(grid,u,X);}
    T2 Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {return Interpolate_Face_V(grid,u,X,axis);}
    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
    {return Extrema_Face_V(grid,u0,u1,X);}
    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
    {return Extrema_Face_V(grid,u0,u1,X,axis);}
    VECTOR<TV2,TV::m> Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {return Gradient_Face_V(grid,u,X);}
    TV2 Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {return Gradient_Face_V(grid,u,X,axis);}
    VECTOR<SM2,TV::m> Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const
    {return Hessian_Face_V(grid,u,X);}
    SM2 Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const
    {return Hessian_Face_V(grid,u,X,axis);}

protected:
    // These are virtual versions of the above; base classes should implement
    // these by calling overridden versions of the routines above.
    virtual TV2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const = 0;
    virtual T2 Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const = 0;
    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 = 0;
    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 = 0;
    virtual VECTOR<TV2,TV::m> Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const = 0;
    virtual TV2 Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const = 0;
    virtual VECTOR<SM2,TV::m> Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const = 0;
    virtual SM2 Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const = 0;

    struct REF
    {
        const T_FACE_ARRAY& u;
        REF(const T_FACE_ARRAY& u,const TV& X):u(u)
        {Set_Reference_Point(u,X);}
        ~REF(){Clear_Reference_Point(u);}
    };
};
//#####################################################################
}
#endif
