//#####################################################################
// Copyright 2002-2004, Ronald Fedkiw.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Advection/ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO.h>
#include <Core/Advection/ADVECTION_SEPARABLE_UNIFORM.h>
namespace PhysBAM{
//#####################################################################
// Function Advection_Solver
//#####################################################################
template<class TV,class T2> void ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<TV,T2>::
Advection_Solver(const int m,const T dx,const ARRAY<T2,IV>& Z,const ARRAY<T,IV>& u,ARRAY<T2,IV>& u_Zx,
        ARRAY<T2>* flux)
{
    assert(!flux);
    T one_over_dx=1/dx,one_over_two_dx=.5*one_over_dx,one_over_three_dx=((T)1/3)*one_over_dx;
    ARRAY<T2,IV> D0(RANGE<IV>(IV(-3),IV(m+3))),D1(RANGE<IV>(IV(-3),IV(m+2))),D2(RANGE<IV>(IV(-3),IV(m+1))); // divided differences
    for(int i=-3;i<m+3;i++) D0(i)=Z(i); 
    if(order>=2) for(int i=-3;i<m+2;i++) D1(i)=(D0(i+1)-D0(i))*one_over_two_dx;     
    if(order==3) for(int i=-3;i<m+1;i++) D2(i)=(D1(i+1)-D1(i))*one_over_three_dx;

    T2 flux_left,flux_right; 
    if(order == 1)
        for(int i=0;i<m;i++)
        {
            if(u(i) > 0)
            {
                flux_left=D0(i-1);flux_right=D0(i);
            }
            else
            {
                flux_left=D0(i);flux_right=D0(i+1);
            }
            u_Zx(i)=u(i)*(flux_right-flux_left)*one_over_dx;
        }
    else if(order == 2)
        for(int i=0;i<m;i++)
        {
            if(u(i) > 0)
            {
                flux_left=ENO(dx,D0(i-1),D1(i-2),D1(i-1));
                flux_right=ENO(dx,D0(i),D1(i-1),D1(i));
            }
            else
            {
                flux_left=ENO(dx,D0(i),-D1(i),-D1(i-1));
                flux_right=ENO(dx,D0(i+1),-D1(i+1),-D1(i));
            }
            u_Zx(i)=u(i)*(flux_right-flux_left)*one_over_dx;
        }
    else if(order == 3)
        for(int i=0;i<m;i++)
        {
            if(u(i) > 0)
            {
                flux_left=ENO(dx,D0(i-1),D1(i-2),D1(i-1),D2(i-3),D2(i-2),D2(i-1));
                flux_right=ENO(dx,D0(i),D1(i-1),D1(i),D2(i-2),D2(i-2),D2(i));
            }
            else
            {
                flux_left=ENO(dx,D0(i),-D1(i),-D1(i-1),D2(i),D2(i-1),D2(i-2));
                flux_right=ENO(dx,D0(i+1),-D1(i+1),-D1(i),D2(i+1),D2(i),D2(i-1));
            }
            u_Zx(i)=u(i)*(flux_right-flux_left)*one_over_dx;
        }
}
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<VECTOR<double,1>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<VECTOR<double,2>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<VECTOR<double,3>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<VECTOR<float,2>,float>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_ENO<VECTOR<float,3>,float>;
} 

