//#####################################################################
// 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_WENO.h>
#include <Core/Advection/ADVECTION_SEPARABLE_UNIFORM.h>
namespace PhysBAM{
//#####################################################################
// Function Advection_Solver
//#####################################################################
template<class TV,class T2> void ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<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;
    ARRAY<T2,IV> D0(RANGE<IV>(IV(-3),IV(m+3))); // 1st divided difference
    for(int i=-3;i<m+3;i++) D0(i)=Z(i);

    if(compute_epsilon)
    {
        epsilon=1e-6*sqr(D0.Max_Abs());
        if(epsilon == 0) epsilon=1e-6;
    } // epsilon=0 implies all v_i=0 and all u_Zx=0

    for(int i=0;i<m;i++)
    {
        if(u(i)>0)
        {
            T2 flux_left=WENO(D0(i-3),D0(i-2),D0(i-1),D0(i),D0(i+1),epsilon),flux_right=WENO(D0(i-2),D0(i-1),D0(i),D0(i+1),D0(i+2),epsilon);
            u_Zx(i)=u(i)*(flux_right-flux_left)*one_over_dx;
        }
        else
        { 
            T2 flux_left=WENO(D0(i+2),D0(i+1),D0(i),D0(i-1),D0(i-2),epsilon),flux_right=WENO(D0(i+3),D0(i+2),D0(i+1),D0(i),D0(i-1),epsilon);
            u_Zx(i)=u(i)*(flux_right-flux_left)*one_over_dx;
        }
    }
}
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<VECTOR<double,1>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<VECTOR<double,2>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<VECTOR<double,3>,double>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<VECTOR<float,2>,float>;
template class ADVECTION_NONCONSERVATIVE_FLUX_FORM_WENO<VECTOR<float,3>,float>;
}

