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

#include <Core/Interpolation/GENERAL_INTERPOLATION_CELL.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_FACE.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_NODE.h>
#include <Core/Interpolation/INTERPOLATION_UNIFORM_FORWARD.h>
namespace PhysBAM{
struct CUBIC_INTERPOLATION_STENCIL
{
    constexpr static int stencil_size=4;
    template<int a,int d,class T> static VECTOR<VECTOR<T,4>,1> Value(T t)
    {
        T s=1-t,u=2-t,r=1+t;
        return {{-(T)1/6*s*t*u,(T).5*r*s*u,(T).5*r*t*u,-(T)1/6*r*s*t}};
    }
    template<int a,int d,class T> static VECTOR<VECTOR<T,4>,1> Diff(T t)
    {
        T s=1-t;
        return {{(T)1/6-(T).5*s*s,(T).5*s*(1-3*t)-1,(T).5*t*(2-3*t)+1,-(T)1/6+(T).5*t*t}};
    }
    template<int a,int d,class T> static VECTOR<VECTOR<T,4>,1> Hess(T t)
    {
        T s=1-t;
        return {{s,3*t-2,1-3*t,t}};
    }
};
template<class TV,class T2>
using CUBIC_INTERPOLATION_NODE=GENERAL_INTERPOLATION_NODE<TV,T2,CUBIC_INTERPOLATION_STENCIL>;

template<class TV,class T2>
using CUBIC_INTERPOLATION_CELL=GENERAL_INTERPOLATION_CELL<TV,T2,CUBIC_INTERPOLATION_STENCIL>;

template<class TV,class T2,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
using CUBIC_INTERPOLATION_FACE=GENERAL_INTERPOLATION_FACE<TV,T2,CUBIC_INTERPOLATION_STENCIL,false,T_FACE_ARRAY>;

template<class TV,class T2,class T_FACE_ARRAY>
struct CUBIC_INTERPOLATION_ALL:CUBIC_INTERPOLATION_NODE<TV,T2>,CUBIC_INTERPOLATION_CELL<TV,T2>,CUBIC_INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY>
{
};
}
#endif
