//#####################################################################
// 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_DEFINITION__
#define __GENERAL_INTERPOLATION_FACE_DEFINITION__

#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_FACE.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_HELPERS.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>
#include <Core/Interpolation/INTERPOLATION_HELPERS.h>
#include <Core/Log/DEBUG_UTILITIES.h>
namespace PhysBAM{

//#####################################################################
// Function Interpolate_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> TV2
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X);
    TV2 r;
    for(int i=0;i<TV::m;i++) r(i)=Interpolate<-1>(stencil,grid,u.Component(i),p(i));
    return r;
}
//#####################################################################
// Function Interpolate_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Interpolate_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> T2
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X,axis);
    return Interpolate<-1>(stencil,grid,u.Component(axis),p);
}
//#####################################################################
// Function Extrema_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X) const -> VECTOR<VECTOR<T2,2>,TV::m>
{
    REF ref0(u0,X),ref1(u1,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u0.Domain_Indices(),X);
    VECTOR<VECTOR<T2,2>,TV::m> r;
    for(int i=0;i<TV::m;i++) r(i)=Extrema<ss,-1>(grid,u0.Component(i),u1.Component(i),p(i));
    return r;
}
//#####################################################################
// Function Extrema_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Extrema_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X,int axis) const -> VECTOR<T2,2>
{
    REF ref0(u0,X),ref1(u1,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u0.Domain_Indices(),X,axis);
    return Extrema<ss,-1>(grid,u0.Component(axis),u1.Component(axis),p);
}
//#####################################################################
// Function Gradient_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> VECTOR<TV2,TV::m>
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X);
    VECTOR<TV2,TV::m> r;
    for(int i=0;i<TV::m;i++) r(i)=Gradient<-1>(stencil,grid,u.Component(i),p(i));
    return r;
}
//#####################################################################
// Function Gradient_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Gradient_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> TV2
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X,axis);
    return Gradient<-1>(stencil,grid,u.Component(axis),p);
}
//#####################################################################
// Function Hessian_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> VECTOR<SM2,TV::m>
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X);
    VECTOR<SM2,TV::m> r;
    for(int i=0;i<TV::m;i++) r(i)=Hessian<-1>(stencil,grid,u.Component(i),p(i));
    return r;
}
//#####################################################################
// Function Hessian_Face
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Hessian_Face(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> SM2
{
    REF ref(u,X);
    auto p=INTERPOLATION_HELPERS::Clamped_Index_Face<ss,const_extrap>(grid,u.Domain_Indices(),X,axis);
    return Hessian<-1>(stencil,grid,u.Component(axis),p);
}
//#####################################################################
// Function Interpolate_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> TV2
{
    return Interpolate_Face(grid,u,X);
}
//#####################################################################
// Function Interpolate_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Interpolate_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> T2
{
    return Interpolate_Face(grid,u,X,axis);
}
//#####################################################################
// Function Extrema_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X) const -> VECTOR<VECTOR<T2,2>,TV::m>
{
    return Extrema_Face(grid,u0,u1,X);
}
//#####################################################################
// Function Extrema_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Extrema_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u0,const T_FACE_ARRAY& u1,const TV& X,int axis) const -> VECTOR<T2,2>
{
    return Extrema_Face(grid,u0,u1,X,axis);
}
//#####################################################################
// Function Gradient_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> VECTOR<TV2,TV::m>
{
    return Gradient_Face(grid,u,X);
}
//#####################################################################
// Function Gradient_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Gradient_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> TV2
{
    return Gradient_Face(grid,u,X,axis);
}
//#####################################################################
// Function Hessian_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X) const -> VECTOR<SM2,TV::m>
{
    return Hessian_Face(grid,u,X);
}
//#####################################################################
// Function Hessian_Face_V
//#####################################################################
template<class TV,class T2,class STENCIL,bool const_extrap,class T_FACE_ARRAY> auto
GENERAL_INTERPOLATION_FACE<TV,T2,STENCIL,const_extrap,T_FACE_ARRAY>::
Hessian_Face_V(const GRID<TV>& grid,const T_FACE_ARRAY& u,const TV& X,int axis) const -> SM2
{
    return Hessian_Face(grid,u,X,axis);
}
}
#endif
