//#####################################################################
// Copyright 2002-2007, Ronald Fedkiw, Geoffrey Irving, Frank Losasso, Avi Robinson-Mosher, Andrew Selle, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class ADVECTION_SEPARABLE_UNIFORM
//#####################################################################
#include <Core/Advection/ADVECTION_SEPARABLE_UNIFORM.h>
#include <Core/Arrays/ARRAYS_ND.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Grids/CELL_ITERATOR.h>
#include <Core/Grids/FACE_ITERATOR.h>
using namespace PhysBAM;
namespace{
// The following construct is necessary over straight function overloading to
// keep MSVC9 from getting confused :(
template<int d> struct UPDATE_ADVECTION_EQUATION_HELPER;
//#####################################################################
// Struct UPDATE_ADVECTION_EQUATION_HELPER<1>
//#####################################################################
template<> struct UPDATE_ADVECTION_EQUATION_HELPER<1>
{
    template<class T_ADVECTION_SEPARABLE_UNIFORM,class T,class T2>
        static void Apply(T_ADVECTION_SEPARABLE_UNIFORM& advection,const GRID<VECTOR<T,1>>& grid,ARRAYS_ND_BASE<T2,VECTOR<int,1> >& Z,
            const ARRAYS_ND_BASE<T2,VECTOR<int,1> >& Z_ghost,const ARRAYS_ND_BASE<VECTOR<T,1>,VECTOR<int,1> >& V,const T dt,const T time,
            ARRAY<T2,FACE_INDEX<1> >* flux_input)
    {
        int i;int m=grid.Numbers_Of_Nodes().x;T dx=grid.dX.x;
        ARRAY<T2,VECTOR<int,1> > rhs(grid.Numbers_Of_Nodes());
        ARRAY<T2,VECTOR<int,1> > Z_1d_x(grid.Node_Indices(3));
        ARRAY<T,VECTOR<int,1> > u_1d(V.domain);
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().x+1);
        for(i=-3;i<m+3;i++) Z_1d_x(i)=Z_ghost(i);
        for(i=V.domain.min_corner.x;i<V.domain.max_corner.x;i++) u_1d(i)=V(i).x;
        advection.Advection_Solver(m,dx,Z_1d_x,u_1d,rhs,flux_input?&flux:0);

        for(i=0;i<m;i++) Z(i)-=dt*rhs(i);

        if(flux_input)
            for(int i=0;i<flux.m;i++)
                (*flux_input)({0,{i}})=flux(i);
    }
};
//#####################################################################
// Struct UPDATE_ADVECTION_EQUATION_HELPER<2>
//#####################################################################
template<> struct UPDATE_ADVECTION_EQUATION_HELPER<2>
{
    template<class T_ADVECTION_SEPARABLE_UNIFORM,class T,class T2>
        static void Apply(T_ADVECTION_SEPARABLE_UNIFORM& advection,const GRID<VECTOR<T,2>>& grid,ARRAYS_ND_BASE<T2,VECTOR<int,2> >& Z,
        const ARRAYS_ND_BASE<T2,VECTOR<int,2> >& Z_ghost,const ARRAYS_ND_BASE<VECTOR<T,2>,VECTOR<int,2> >& V,const T dt,const T time,
            ARRAY<T2,FACE_INDEX<2> >* flux_input)
    {
        int i,j;int m=grid.Numbers_Of_Nodes().x,n=grid.Numbers_Of_Nodes().y;T dx=grid.dX.x,dy=grid.dX.y;ARRAY<T2,VECTOR<int,2> > rhs(grid.Numbers_Of_Nodes());

        RANGE<VECTOR<int,1> > R;
        R.min_corner.x=0;
        R.max_corner.x=m;
        {ARRAY<T2,VECTOR<int,1> > Z_1d_x(R.Thickened(3)),u_Zx_1d(R);ARRAY<T,VECTOR<int,1> > u_1d(V.domain.Dimension_Range(0));
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().x+1);
        for(j=0;j<n;j++){
            for(i=-3;i<m+3;i++) Z_1d_x(i)=Z_ghost(i,j);
            for(i=V.domain.min_corner.x;i<V.domain.max_corner.x;i++) u_1d(i)=V(i,j).x;
            advection.Advection_Solver(m,dx,Z_1d_x,u_1d,u_Zx_1d,flux_input?&flux:0);
            for(i=0;i<m;i++) rhs(i,j)=u_Zx_1d(i);

            if(flux_input)
                for(int i=0;i<flux.m;i++)
                    (*flux_input)({0,{i,j}})=flux(i);
        }}

        R.max_corner.x=n;
        {ARRAY<T2,VECTOR<int,1> > Z_1d_y(R.Thickened(3)),v_Zy_1d(R);ARRAY<T,VECTOR<int,1> > v_1d(V.domain.Dimension_Range(1));
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().y+1);
        for(i=0;i<m;i++){
            for(j=-3;j<n+3;j++) Z_1d_y(j)=Z_ghost(i,j);
            for(j=V.domain.min_corner.y;j<V.domain.max_corner.y;j++) v_1d(j)=V(i,j).y;
            advection.Advection_Solver(n,dy,Z_1d_y,v_1d,v_Zy_1d,flux_input?&flux:0);
            for(j=0;j<n;j++) rhs(i,j)+=v_Zy_1d(j);

            if(flux_input)
                for(int j=0;j<flux.m;j++)
                    (*flux_input)({1,{i,j}})=flux(j);
        }}

        for(i=0;i<m;i++) for(j=0;j<n;j++) Z(i,j)-=dt*rhs(i,j);
    }
};
//#####################################################################
// Struct UPDATE_ADVECTION_EQUATION_HELPER<3>
//#####################################################################
template<> struct UPDATE_ADVECTION_EQUATION_HELPER<3>
{
    template<class T_ADVECTION_SEPARABLE_UNIFORM,class T,class T2>
        static void Apply(T_ADVECTION_SEPARABLE_UNIFORM& advection,const GRID<VECTOR<T,3>>& grid,ARRAYS_ND_BASE<T2,VECTOR<int,3> >& Z,
        const ARRAYS_ND_BASE<T2,VECTOR<int,3> >& Z_ghost,const ARRAYS_ND_BASE<VECTOR<T,3>,VECTOR<int,3> >& V,const T dt,const T time,
            ARRAY<T2,FACE_INDEX<3> >* flux_input)
    {
        int i,j,k;
        int m=grid.Numbers_Of_Nodes().x,n=grid.Numbers_Of_Nodes().y,mn=grid.Numbers_Of_Nodes().z;
        T dx=grid.dX.x,dy=grid.dX.y,dz=grid.dX.z;
        ARRAY<T2,VECTOR<int,3> > rhs(grid.Numbers_Of_Nodes());

        RANGE<VECTOR<int,1> > R;
        R.min_corner.x=0;
        R.max_corner.x=m;
        {ARRAY<T2,VECTOR<int,1> > Z_1d_x(R.Thickened(3)),u_Zx_1d(R);ARRAY<T,VECTOR<int,1> > u_1d(V.domain.Dimension_Range(0));
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().x+1);
        for(j=0;j<n;j++) for(k=0;k<mn;k++){
            for(i=-3;i<m+3;i++) Z_1d_x(i)=Z_ghost(i,j,k);
            for(i=V.domain.min_corner.x;i<V.domain.max_corner.x;i++) u_1d(i)=V(i,j,k).x;
            advection.Advection_Solver(m,dx,Z_1d_x,u_1d,u_Zx_1d,flux_input?&flux:0);
            for(i=0;i<m;i++) rhs(i,j,k)=u_Zx_1d(i);

            if(flux_input)
                for(int i=0;i<flux.m;i++)
                    (*flux_input)({0,{i,j,k}})=flux(i);
        }}

        R.max_corner.x=n;
        {ARRAY<T2,VECTOR<int,1> > Z_1d_y(R.Thickened(3)),v_Zy_1d(R);ARRAY<T,VECTOR<int,1> > v_1d(V.domain.Dimension_Range(1));
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().y+1);
        for(i=0;i<m;i++) for(k=0;k<mn;k++){
            for(j=-3;j<n+3;j++) Z_1d_y(j)=Z_ghost(i,j,k);
            for(j=V.domain.min_corner.y;j<V.domain.max_corner.y;j++) v_1d(j)=V(i,j,k).y;
            advection.Advection_Solver(n,dy,Z_1d_y,v_1d,v_Zy_1d,flux_input?&flux:0);
            for(j=0;j<n;j++) rhs(i,j,k)+=v_Zy_1d(j);

            if(flux_input)
                for(int j=0;j<flux.m;j++)
                    (*flux_input)({1,{i,j,k}})=flux(j);
        }}

        R.max_corner.x=mn;
        {ARRAY<T2,VECTOR<int,1> > Z_1d_z(R.Thickened(3)),w_Zz_1d(R);ARRAY<T,VECTOR<int,1> > w_1d(V.domain.Dimension_Range(2));
        ARRAY<T2> flux;
        if(flux_input) flux.Resize(grid.Numbers_Of_Nodes().z+1);
        for(i=0;i<m;i++) for(j=0;j<n;j++){
            for(k=-3;k<mn+3;k++) Z_1d_z(k)=Z_ghost(i,j,k);
            for(k=V.domain.min_corner.z;k<V.domain.max_corner.z;k++) w_1d(k)=V(i,j,k).z;
            advection.Advection_Solver(mn,dz,Z_1d_z,w_1d,w_Zz_1d,flux_input?&flux:0);
            for(k=0;k<mn;k++) rhs(i,j,k)+=w_Zz_1d(k);

            if(flux_input)
                for(int k=0;k<flux.m;k++)
                    (*flux_input)({2,{i,j,k}})=flux(k);
        }}

        for(RANGE_ITERATOR<3> it(grid.Node_Indices());it.Valid();it.Next()) Z(it.index)-=dt*rhs(it.index);
    }
};
}
//#####################################################################
// Function Update_Advection_Equation_Node
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Update_Advection_Equation_Node(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,const ARRAY<TV,IV>& V,const T dt,const T time,
    const ARRAY<T2,IV>* Z_min_ghost,const ARRAY<T2,IV>* Z_max_ghost,ARRAY<T2,IV>* Z_min,ARRAY<T2,IV>* Z_max,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux)
{
    assert(!Z_min && !Z_max);
    UPDATE_ADVECTION_EQUATION_HELPER<TV::m>::Apply(*this,grid,Z,Z_ghost,V,dt,time,flux);
}
//#####################################################################
// Function Update_Advection_Equation_Cell_Lookup
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Update_Advection_Equation_Cell_Lookup(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,const T_FACE_LOOKUP& V,const T dt,const T time,
    const ARRAY<T2,IV>* Z_min_ghost,const ARRAY<T2,IV>* Z_max_ghost,ARRAY<T2,IV>* Z_min,ARRAY<T2,IV>* Z_max,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux)
{
    assert(!Z_min && !Z_max);

    ARRAY<TV,IV> V_cell(grid.Cell_Indices(3));T_AVERAGING averaging;

    for(CELL_ITERATOR<TV> iterator(grid,3);iterator.Valid();iterator.Next()) V_cell(iterator.Cell_Index())=averaging.Face_To_Cell_Vector(grid,iterator.Cell_Index(),V);

    UPDATE_ADVECTION_EQUATION_HELPER<TV::m>::Apply(*this,grid.Get_Regular_Grid_At_MAC_Positions(),Z,Z_ghost,V_cell,dt,time,flux);
}
//#####################################################################
// Function Update_Advection_Equation_Cell
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Update_Advection_Equation_Face_Lookup(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& Z,const T_FACE_LOOKUP& Z_ghost,const T_FACE_LOOKUP& V,
    const T dt,const T time,const T_FACE_LOOKUP* Z_min_ghost,const T_FACE_LOOKUP* Z_max_ghost,ARRAY<T,FACE_INDEX<TV::m> >* Z_min,ARRAY<T,FACE_INDEX<TV::m> >* Z_max,
        VECTOR<ARRAY<T2,FACE_INDEX<TV::m> >,TV::m>* flux)
{
    assert(!Z_min && !Z_max);

    ARRAY<TV,FACE_INDEX<TV::m> > V_face(grid.Cell_Indices(3));
    T_AVERAGING averaging;

    for(FACE_ITERATOR<TV> iterator(grid,3);iterator.Valid();iterator.Next()) V_face(iterator.Full_Index())=averaging.Face_To_Face_Vector(grid,iterator.Full_Index(),V);

    for(int i=0;i<TV::m;i++){
        GRID<TV> node_grid(grid.Get_Face_Grid(i));
        UPDATE_ADVECTION_EQUATION_HELPER<TV::m>::Apply(*this,node_grid,Z.Component(i),Z_ghost.Component(i),V_face.Component(i),dt,time,flux?&(*flux)(i):0);}
}
//#####################################################################
// Function Set_Order
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Set_Order(const int order_input)
{
    PHYSBAM_FUNCTION_IS_NOT_DEFINED();
}
//#####################################################################
// Function Set_Epsilon
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Set_Epsilon(const T epsilon_input)
{
    PHYSBAM_FUNCTION_IS_NOT_DEFINED();
}
//#####################################################################
// Function Set_Axis
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Set_Axis(const int axis)
{
    PHYSBAM_FUNCTION_IS_NOT_DEFINED();
}
//#####################################################################
// Function Compute_Epsilon
//#####################################################################
template<class TV,class T2,class T_AVERAGING> void ADVECTION_SEPARABLE_UNIFORM<TV,T2,T_AVERAGING>::
Compute_Epsilon()
{
    PHYSBAM_FUNCTION_IS_NOT_DEFINED();
}
//#####################################################################
namespace PhysBAM{
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<float,1>,float>;
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<float,2>,float>;
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<float,3>,float>;
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<double,1>,double>;
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<double,2>,double>;
template class ADVECTION_SEPARABLE_UNIFORM<VECTOR<double,3>,double>;
}
