//#####################################################################
// Copyright 2022, Ritoban Roy-Chowdhury
// 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_FLUXES.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/LINEAR_INTERPOLATION.h>
namespace PhysBAM{

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