//##################################################################### 
// 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_INTERP.h>
#include <Core/Interpolation/DIV_FREE_INTERPOLATION.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_AXIS_FACE_DEFINITION.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_HELPERS.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>
#include <Core/Interpolation/INTERPOLATION_HELPERS.h>
#include <Core/Interpolation/INTERPOLATION_SPLINES.h>
#include <Core/Symbolics/SYMBOLIC_NUMBER.h>
#include <Core/Vectors/VECTOR.h>

namespace PhysBAM {

// interpolating polynomial
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<2>::
Value(T t) -> VECTOR<VECTOR<T,2+(a==d)>,3> 
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline(IS::B2,u);
        auto p1=Evaluate_Spline(IS::C3,u);
        auto p2=Evaluate_Spline(IS::D4,u);
        return {p0,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline(IS::B1,u);
        auto p1=Evaluate_Spline(IS::C2,u);
        auto p2=Evaluate_Spline(IS::D3,u);
        return {p0,p2,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<2>::
Diff(T t) -> VECTOR<VECTOR<T,2+(a==d)>,3>
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Diff(IS::B2,u);
        auto p1=Evaluate_Spline_Diff(IS::C3,u);
        auto p2=Evaluate_Spline_Diff(IS::D4,u);
        return {p0,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<2>(t);
        auto p0=Evaluate_Spline_Diff(IS::B1,u);
        auto p1=Evaluate_Spline_Diff(IS::C2,u);
        auto p2=Evaluate_Spline_Diff(IS::D3,u);
        return {p0,p2,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<2>::
Hess(T t) -> VECTOR<VECTOR<T,2+(a==d)>,3>
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<2>(t);
        auto p0=Evaluate_Spline_Hess(IS::B2,u);
        auto p1=Evaluate_Spline_Hess(IS::C3,u);
        auto p2=Evaluate_Spline_Hess(IS::D4,u);
        return {p0,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<1>(t);
        auto p0=Evaluate_Spline_Hess(IS::B1,u);
        auto p1=Evaluate_Spline_Hess(IS::C2,u);
        auto p2=Evaluate_Spline_Hess(IS::D3,u);
        return {p0,p2,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<3>::
Value(T t) -> VECTOR<VECTOR<T,2+(a==d)>,4>
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<4>(t);
        auto p0=Evaluate_Spline(IS::B2,u);
        auto p1=Evaluate_Spline(IS::C3,u);
        auto p2=Evaluate_Spline(IS::D4,u);
        return {p0,-4*p1,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline(IS::B1,u);
        auto p1=Evaluate_Spline(IS::C2,u);
        auto p2=Evaluate_Spline(IS::D3,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p2,p1};
        else return {p0,p2,p1,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<3>::
Diff(T t) -> VECTOR<VECTOR<T,2+(a==d)>,4>
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<3>(t);
        auto p0=Evaluate_Spline_Diff(IS::B2,u);
        auto p1=Evaluate_Spline_Diff(IS::C3,u);
        auto p2=Evaluate_Spline_Diff(IS::D4,u);
        return {p0,-4*p1,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<2>(t);
        auto p0=Evaluate_Spline_Diff(IS::B1,u);
        auto p1=Evaluate_Spline_Diff(IS::C2,u);
        auto p2=Evaluate_Spline_Diff(IS::D3,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p2,p1};
        else return {p0,p2,p1,p1};
    }
}
template<int a,int d,class T> auto
DIV_FREE_INTERPOLATION_INTERP_STENCIL<3>::
Hess(T t) -> VECTOR<VECTOR<T,2+(a==d)>,4>
{
    if constexpr(a==d)
    {
        auto u=Make_Powers<2>(t);
        auto p0=Evaluate_Spline_Hess(IS::B2,u);
        auto p1=Evaluate_Spline_Hess(IS::C3,u);
        auto p2=Evaluate_Spline_Hess(IS::D4,u);
        return {p0,-4*p1,-4*p1,-4*p2};
    }
    else
    {
        auto u=Make_Powers<1>(t);
        auto p0=Evaluate_Spline_Hess(IS::B1,u);
        auto p1=Evaluate_Spline_Hess(IS::C2,u);
        auto p2=Evaluate_Spline_Hess(IS::D3,u);
        if constexpr((a-d-1)%3==0) return {p0,p1,p2,p1};
        else return {p0,p2,p1,p1};
    }
}
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,2>,double,DIV_FREE_INTERPOLATION_INTERP_STENCIL<2>,false,ARRAY<double,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,3>,double,DIV_FREE_INTERPOLATION_INTERP_STENCIL<3>,false,ARRAY<double,FACE_INDEX<3> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,2>,float,DIV_FREE_INTERPOLATION_INTERP_STENCIL<2>,false,ARRAY<float,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,3>,float,DIV_FREE_INTERPOLATION_INTERP_STENCIL<3>,false,ARRAY<float,FACE_INDEX<3> > >;
}
