//#####################################################################
// Copyright 2002-2010, Ronald Fedkiw, Geoffrey Irving, Nipun Kwatra, Michael Lentine, Frank Losasso, Andrew Selle.
// 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_SEMI_LAGRANGIAN_UNIFORM_BETA.h>
#include <Core/Grids/CELL_ITERATOR.h>
#include <Core/Grids/FACE_ITERATOR.h>
#include <Core/Grids/GRID.h>
#include <Core/Grids/NODE_ITERATOR.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class TV,class T2,class T_AVERAGING,class T_INTERPOLATION> ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<TV,T2,T_AVERAGING,T_INTERPOLATION>::
ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA()
{}
//#####################################################################
// Function Update_Advection_Equation_Node
//#####################################################################
template<class TV,class T2,class T_AVERAGING,class T_INTERPOLATION> void ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<TV,T2,T_AVERAGING,T_INTERPOLATION>::
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(!flux);
    RANGE<IV> domain=grid.Node_Indices();domain.max_corner+=IV::All_Ones_Vector();
    T_INTERPOLATION interpolation;
    if(Z_min && Z_max) for(NODE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV node=iterator.Node_Index();TV X=iterator.Location()-dt*V(node);
        Z(node)=interpolation.Interpolate_Node(grid,Z_ghost,X);
        VECTOR<T2,2> extrema=interpolation.Extrema_Node(grid,*Z_min_ghost,*Z_max_ghost,X);
        (*Z_min)(node)=extrema.x;(*Z_max)(node)=extrema.y;}
    else for(NODE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){IV node=iterator.Node_Index();
        Z(node)=interpolation.Interpolate_Node(grid,Z_ghost,iterator.Location()-dt*V(node));}
}
//#####################################################################
// Function Update_Advection_Equation_Cell_Lookup
//#####################################################################
template<class TV,class T2,class T_AVERAGING,class T_INTERPOLATION> void ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<TV,T2,T_AVERAGING,T_INTERPOLATION>::
Update_Advection_Equation_Cell_Lookup(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,
        const T_FACE_LOOKUP& face_velocities,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(!flux);
    T_INTERPOLATION interpolation;
    T_AVERAGING averaging;
    if(Z_min && Z_max) for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV cell=iterator.Cell_Index();
        TV X=iterator.Location()-dt*averaging.Face_To_Cell_Vector(grid,cell,face_velocities);
        Z(cell)=interpolation.Interpolate_Cell(grid,Z_ghost,X);
        VECTOR<T2,2> extrema=interpolation.Extrema_Cell(grid,*Z_min_ghost,*Z_max_ghost,X);
        (*Z_min)(cell)=extrema.x;(*Z_max)(cell)=extrema.y;}
    else for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV cell=iterator.Cell_Index();
        TV cell_velocity(averaging.Face_To_Cell_Vector(grid,cell,face_velocities));
        cell_velocity*=dt; // TODO: gcc 4.1.2 compiler bug workaround
        Z(cell)=interpolation.Interpolate_Cell(grid,Z_ghost,iterator.Location()-cell_velocity);}
}
//#####################################################################
// Function Update_Advection_Equation_Face_Lookup
//#####################################################################
template<class TV,class T2,class T_AVERAGING,class T_INTERPOLATION> void ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<TV,T2,T_AVERAGING,T_INTERPOLATION>::
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& face_velocities,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(!flux);
    if constexpr(is_same<typename SCALAR_POLICY<T2>::TYPE,T2>::value)
    {
        T_INTERPOLATION interpolation;
        T_AVERAGING averaging;
        if(Z_min && Z_max)
        {
            for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
                IV face=iterator.Face_Index();int axis=iterator.Axis();TV X=iterator.Location()-dt*averaging.Face_To_Face_Vector(grid,axis,face,face_velocities);
                Z.Component(axis)(face)=interpolation.Interpolate_Face(grid,Starting_Point_Face(Z_ghost,axis,face),X,axis);
                VECTOR<T,2> extrema=interpolation.Extrema_Face(grid,Starting_Point_Face(*Z_min_ghost,axis,face),Starting_Point_Face(*Z_max_ghost,axis,face),X,axis);
                (*Z_min).Component(axis)(face)=extrema.x;(*Z_max).Component(axis)(face)=extrema.y;}
        }
        else
            for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){IV face=iterator.Face_Index();int axis=iterator.Axis();
                Z.Component(axis)(face)=interpolation.Interpolate_Face(grid,Starting_Point_Face(Z_ghost,axis,face),
                    iterator.Location()-dt*averaging.Face_To_Face_Vector(grid,axis,face,face_velocities),axis);}
    }
    else PHYSBAM_FATAL_ERROR();
}
