//#####################################################################
// Copyright 2003-2021, Ronald Fedkiw, Geoffrey Irving, Nipun Kwatra, Duc Nguyen, Craig Schroeder, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class CUBIC_MN_INTERPOLATION 
//#####################################################################
#ifndef __CUBIC_MN_INTERPOLATION__
#define __CUBIC_MN_INTERPOLATION__

#include <Core/Arrays/ARRAY.h>
#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{

template<class T>
class CUBIC_MN_INTERPOLATION
{
private:
    T m0,m1,m3,n0,n1,n2,n3; // used for acceleration

    template<class T2,int d>
    static void sign_check(VECTOR<T2,d>& di,const VECTOR<T2,d>& dq)
    {
        for(int i=0;i<d;i++) if(di(i)*dq(i)<0) di(i)=0;
    }
    template<class T2,int d>
    static void sign_check(MATRIX<T2,d>& di,const MATRIX<T2,d>& dq)
    {
        for(int i=0;i<d;i++) for(int j=0;j<d;j++) if(di(i,j)*dq(i,j)<0) di(i,j)=0;
    }
    template<class T2,int d>
    static void sign_check(SYMMETRIC_MATRIX<T2,d>& di,const SYMMETRIC_MATRIX<T2,d>& dq)
    {
        for(int i=0;i<d;i++)
            for(int j=0;j<=i;j++)
                if(di.Element_Lower(i,j)*dq.Element_Lower(i,j)<0)
                    di.Element_Lower(i,j)=0;
    }
    template<class T2>
    static void sign_check(T2& di,const T2& dq)
    {
        if(di*dq<0) di=0;
    }
public:

    CUBIC_MN_INTERPOLATION();
    ~CUBIC_MN_INTERPOLATION();
    void Set_Parameters(const T b=1./3,const T c=1./3);
    ARRAY<T> Cubic_MN_Weights(const T alpha) const;

    template<class T2> T2
    Cubic_MN(const T2& u_0,const T2& u_1,const T2& u_2,const T2& u_3,const T alpha,const bool monotonic) const
    {
        if(monotonic) return Cubic_MN_Monotonic(u_0,u_1,u_2,u_3,alpha);
        T temp=abs(1+alpha);
        T2 result=u_0*(((n0*temp+n1)*temp+n2)*temp+n3);
        temp=abs(alpha);result+=u_1*((m0*temp+m1)*temp*temp+m3);
        temp=abs(1-alpha);result+=u_2*((m0*temp+m1)*temp*temp+m3);
        temp=abs(2-alpha);result+=u_3*(((n0*temp+n1)*temp+n2)*temp+n3);
        return result;
    }

    template<class T2> T2
    Cubic_MN_Monotonic(const T2& u_0,const T2& u_1,const T2& u_2,const T2& u_3,const T alpha) const
    {
        T x2=alpha*alpha;
        T x3=x2*alpha;
        T2 dq=u_2-u_1;
        T2 di=(T).5*(u_2-u_0);
        sign_check(di,dq);
        T2 dip1=(T).5*(u_3-u_1);
        sign_check(dip1,dq);
        return u_1+di*alpha+((T)3.*dq-(T)2.*di-dip1)*x2+(-(T)2.*dq+di+dip1)*x3;
    }
    
//#####################################################################
};

template<class T>
struct CUBIC_MN_INTERPOLATION_STENCIL
{
    CUBIC_MN_INTERPOLATION<T> cub;
    bool monotonic=false;
    constexpr static int stencil_size=4;
    template<class T2> VECTOR<T2,1> Value(T t,const T2* x) const
    {return {cub.Cubic_MN(x[0],x[1],x[2],x[3],t,monotonic)};}
    template<class T2> VECTOR<T2,1> Diff(T t,const T2* x) const
    {PHYSBAM_FATAL_ERROR();}
    template<class T2> VECTOR<T2,1> Hess(T t,const T2* x) const
    {PHYSBAM_FATAL_ERROR();}
};

template<class TV,class T2>
using CUBIC_MN_INTERPOLATION_NODE=GENERAL_INTERPOLATION_NODE<TV,T2,CUBIC_MN_INTERPOLATION_STENCIL<typename TV::SCALAR> >;

template<class TV,class T2>
using CUBIC_MN_INTERPOLATION_CELL=GENERAL_INTERPOLATION_CELL<TV,T2,CUBIC_MN_INTERPOLATION_STENCIL<typename TV::SCALAR> >;

template<class TV,class T2,class T_FACE_ARRAY=ARRAY<T2,FACE_INDEX<TV::m> > >
using CUBIC_MN_INTERPOLATION_FACE=GENERAL_INTERPOLATION_FACE<TV,T2,CUBIC_MN_INTERPOLATION_STENCIL<typename TV::SCALAR>,false,T_FACE_ARRAY>;

template<class TV,class T2,class T_FACE_ARRAY>
struct CUBIC_MN_INTERPOLATION_ALL:CUBIC_MN_INTERPOLATION_NODE<TV,T2>,CUBIC_MN_INTERPOLATION_CELL<TV,T2>,CUBIC_MN_INTERPOLATION_FACE<TV,T2,T_FACE_ARRAY>
{
    typedef CUBIC_MN_INTERPOLATION_NODE<TV,T2> N;
    typedef CUBIC_MN_INTERPOLATION_CELL<TV,T2> C;
    typedef CUBIC_MN_INTERPOLATION_FACE<TV,T2> F;
    typedef typename TV::SCALAR T;
    
    void Set_Sharpness(T sharpness_input=2./3)
    {
        N::stencil.Set_Sharpness(sharpness_input);
        C::stencil.Set_Sharpness(sharpness_input);
        F::stencil.Set_Sharpness(sharpness_input);
    }

    void Set_Parameters(const T b_input=1./3,const T c_input=1./3)
    {
        N::stencil.Set_Paremeters(b_input,c_input);
        C::stencil.Set_Paremeters(b_input,c_input);
        F::stencil.Set_Paremeters(b_input,c_input);
    }

    void Set_Monotonic(bool monototic=false)
    {
        N::stencil.Set_Monotonic(monototic);
        C::stencil.Set_Monotonic(monototic);
        F::stencil.Set_Monotonic(monototic);
    }
};

}
#endif
