//#####################################################################
// 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.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>
namespace PhysBAM{

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