//#####################################################################
// Copyright 2022, Ritoban Roy-Chowdhury, 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/CURL_FREE_INTERPOLATION_C1.h>
#include <Core/Arrays/ARRAY.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Interpolation/GENERAL_INTERPOLATION_AXIS_FACE_DEFINITION.h>
#include <Core/Interpolation/INTERPOLATION_FACE.h>
#include <Core/Interpolation/INTERPOLATION_HELPERS.h>
#include <Core/Interpolation/INTERPOLATION_SPLINES.h>
namespace PhysBAM{

template<int a,int d,class T>
VECTOR<VECTOR<T,4+(a==d)>,1> CURL_FREE_INTERPOLATION_C1_STENCIL::
Value(T t)
{
    T s=1-t,s2=s*s,t2=t*t;
    if constexpr(a==d) return {{0,(T).5*s2,s*t+(T).5,(T).5*t2,0}};
    else return {{(T)1/6*s2*s,(T).5*t2*(t-2)+(T)2/3,-(T).5*s2*(t+1)+(T)2/3,(T)1/6*t2*t}};
}
template<int a,int d,class T>
VECTOR<VECTOR<T,4+(a==d)>,1> CURL_FREE_INTERPOLATION_C1_STENCIL::
Diff(T t)
{
    if constexpr(a==d) return {{0,t-1,1-2*t,t,0}};
    else
    {
        T s=1-t;
        return {{-(T).5*s*s,(T).5*t*(3*t-4),(T).5*s*(3*t+1),(T).5*t*t}};
    }
}
template<int a,int d,class T>
VECTOR<VECTOR<T,4+(a==d)>,1> CURL_FREE_INTERPOLATION_C1_STENCIL::
Hess(T t)
{
    if constexpr(a==d) return {{0,1,-2,1,0}};
    else return {{1-t,3*t-2,1-3*t,t}};
}
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,2>,double,CURL_FREE_INTERPOLATION_C1_STENCIL,false,ARRAY<double,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<double,3>,double,CURL_FREE_INTERPOLATION_C1_STENCIL,false,ARRAY<double,FACE_INDEX<3> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,2>,float,CURL_FREE_INTERPOLATION_C1_STENCIL,false,ARRAY<float,FACE_INDEX<2> > >;
template class GENERAL_INTERPOLATION_AXIS_FACE<VECTOR<float,3>,float,CURL_FREE_INTERPOLATION_C1_STENCIL,false,ARRAY<float,FACE_INDEX<3> > >;
}
