//#####################################################################
// 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.
//#####################################################################
#include <Core/Interpolation/DIV_FREE_INTERPOLATION_4TH_ORDER.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_AXIS_FACE_DEFINITION.h>
#include <Core/Interpolation/INTERPOLATION_SPLINES.h>

namespace PhysBAM{

template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<2>::
Value(T t) -> VECTOR<VECTOR<T,6+(a==d)>,2>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline(IS4::B2p,u);
        auto p1=Evaluate_Spline(IS4::I4,u);
        return {p1*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline(IS::B1pp,u);
        auto p1=Evaluate_Spline(IS::I3p,u);
        return {p0,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<2>::
Diff(T t) -> VECTOR<VECTOR<T,6+(a==d)>,2>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline_Diff(IS4::B2p,u);
        auto p1=Evaluate_Spline_Diff(IS4::I4,u);
        return {p1*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Diff(IS::B1pp,u);
        auto p1=Evaluate_Spline_Diff(IS::I3p,u);
        return {p0,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<2>::
Hess(T t) -> VECTOR<VECTOR<T,6+(a==d)>,2>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline_Hess(IS4::B2p,u);
        auto p1=Evaluate_Spline_Hess(IS4::I4,u);
        return {p1*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Hess(IS::B1pp,u);
        auto p1=Evaluate_Spline_Hess(IS::I3p,u);
        return {p0,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<3>::
Value(T t) -> VECTOR<VECTOR<T,6+(a==d)>,3>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline(IS4::B2p,u);
        auto p1=Evaluate_Spline(IS4::J4,u);
        return {p1*k0,p0*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline(IS::B1pp,u);
        auto p1=Evaluate_Spline(IS::J3p,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p0};
        else return {p0,p0,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<3>::
Diff(T t) -> VECTOR<VECTOR<T,6+(a==d)>,3>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline_Diff(IS4::B2p,u);
        auto p1=Evaluate_Spline_Diff(IS4::J4,u);
        return {p1*k0,p0*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Diff(IS::B1pp,u);
        auto p1=Evaluate_Spline_Diff(IS::J3p,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p0};
        else return {p0,p0,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<3>::
Hess(T t) -> VECTOR<VECTOR<T,6+(a==d)>,3>
{
    if constexpr(a==d)
    {
        T k0=1./6;
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline_Hess(IS4::B2p,u);
        auto p1=Evaluate_Spline_Hess(IS4::J4,u);
        return {p1*k0,p0*k0,p0*k0};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Hess(IS::B1pp,u);
        auto p1=Evaluate_Spline_Hess(IS::J3p,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p0};
        else return {p0,p0,p1};
    }
}
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,2>,double,DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<2>,false,ARRAY<double,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,3>,double,DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<3>,false,ARRAY<double,FACE_INDEX<3> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,2>,float,DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<2>,false,ARRAY<float,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,3>,float,DIV_FREE_INTERPOLATION_4TH_ORDER_STENCIL<3>,false,ARRAY<float,FACE_INDEX<3> > >;
}
