//#####################################################################
// Copyright 2006, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class ADVECTION_MACCORMACK_UNIFORM
//#####################################################################
#include <Core/Advection/ADVECTION_MACCORMACK_UNIFORM.h>
#include <Core/Advection/ADVECTION_SEMI_LAGRANGIAN_UNIFORM.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Boundaries/BOUNDARY.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/Interpolation/FACE_LOOKUP_UNIFORM.h>
namespace PhysBAM{
//#####################################################################
// Constructor
//#####################################################################
template<class TV,class T2,class T_NESTED_ADVECTION> ADVECTION_MACCORMACK_UNIFORM<TV,T2,T_NESTED_ADVECTION>::
ADVECTION_MACCORMACK_UNIFORM(T_NESTED_ADVECTION& nested_advection_input,const BOUNDARY<TV,T>& boundary,const ARRAY<bool,IV>* node_mask_input,const ARRAY<bool,IV>* cell_mask_input,const ARRAY<bool,FACE_INDEX<TV::m> >* face_mask_input)
    :node_mask(node_mask_input),cell_mask(cell_mask_input),
    face_mask(face_mask_input),nested_advection(nested_advection_input),
    clamp_extrema(true),ensure_second_order(false),boundary(&boundary)
{
}
//#####################################################################
// Function Update_Advection_Equation_Node
//#####################################################################
template<class TV,class T2,class T_NESTED_ADVECTION> void ADVECTION_MACCORMACK_UNIFORM<TV,T2,T_NESTED_ADVECTION>::
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_input,const ARRAY<T2,IV>* Z_max_ghost_input,ARRAY<T2,IV>* Z_min_input,ARRAY<T2,IV>* Z_max_input,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux)
{
    assert(!flux);
    assert(node_mask && !Z_min_input && !Z_max_input);
    ARRAY<TV,IV> negative_V(grid.Node_Indices(0),no_init);
    negative_V.Copy((T)-1,V);
    ARRAY<T2,IV> Z_temp=Z,Z_min_ghost=Z_ghost,Z_max_ghost=Z_ghost,Z_min(grid.Node_Indices(),no_init),Z_max(grid.Node_Indices(),no_init);
    int number_of_ghost_cells=Z.Domain_Indices().Minimum_Corner()(0)-Z_ghost.Domain_Indices().Minimum_Corner()(0);
    nested_advection.Update_Advection_Equation_Node(grid,Z_temp,Z_ghost,V,dt,time,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp now has time n+1 data
    ARRAY<T2,IV> Z_forward_ghost(grid.Node_Indices(number_of_ghost_cells),no_init);
    boundary->Fill_Ghost_Nodes(grid,Z_temp,Z_forward_ghost,number_of_ghost_cells);
    boundary->Fill_Ghost_Nodes(grid,Z_min,Z_min_ghost,number_of_ghost_cells);
    boundary->Fill_Ghost_Nodes(grid,Z_max,Z_max_ghost,number_of_ghost_cells);
    nested_advection.Update_Advection_Equation_Node(grid,Z_temp,Z_forward_ghost,negative_V,dt,time+dt,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp has time n data
    RANGE<IV> domain=grid.Node_Indices();domain.max_corner+=IV::All_Ones_Vector();
    for(NODE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV index=iterator.Node_Index();
        if((*node_mask)(index)){
            T2 Z_local=Z_forward_ghost(index)+(T).5*(Z(index)-Z_temp(index));
            if(clamp_extrema) Z_local=clamp(Z_local,Z_min(index),Z_max(index));
            Z(index)=Z_local;}
        else Z(index)=Z_forward_ghost(index);}
}
//#####################################################################
// Function Update_Advection_Equation_Cell_Lookup
//#####################################################################
template<class TV,class T2,class T_NESTED_ADVECTION> void ADVECTION_MACCORMACK_UNIFORM<TV,T2,T_NESTED_ADVECTION>::
Update_Advection_Equation_Cell_Lookup(const GRID<TV>& grid,ARRAY<T2,IV>& Z,const ARRAY<T2,IV>& Z_ghost,const ARRAY<typename TV::SCALAR,FACE_INDEX<TV::m> >& face_velocities,
    const T dt,const T time,const ARRAY<T2,IV>* Z_min_ghost_input,const ARRAY<T2,IV>* Z_max_ghost_input,ARRAY<T2,IV>* Z_min_input,ARRAY<T2,IV>* Z_max_input,
        ARRAY<T2,FACE_INDEX<TV::m> >* flux)
{
    assert(!flux);
    assert(cell_mask && !Z_min_input && !Z_max_input);

    ARRAY<T,FACE_INDEX<TV::m> > negative_V(face_velocities.Domain_Indices());
    negative_V.Copy(-1,face_velocities);

    ARRAY<T2,IV> Z_temp=Z,
        Z_min_ghost=Z_ghost,
        Z_max_ghost=Z_ghost,
        Z_min(grid.Cell_Indices(),no_init),
        Z_max(grid.Cell_Indices(),no_init);

    int number_of_ghost_cells=Z.Domain_Indices().Minimum_Corner()(0)-Z_ghost.Domain_Indices().Minimum_Corner()(0);
    
    nested_advection.Update_Advection_Equation_Cell(grid,Z_temp,Z_ghost,face_velocities,dt,time,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp now has time n+1 data
    
    ARRAY<T2,IV> Z_forward_ghost(grid.Cell_Indices(number_of_ghost_cells),no_init);
    boundary->Fill_Ghost_Cells(grid,Z_temp,Z_forward_ghost,number_of_ghost_cells);
    boundary->Fill_Ghost_Cells(grid,Z_min,Z_min_ghost,number_of_ghost_cells);
    boundary->Fill_Ghost_Cells(grid,Z_max,Z_max_ghost,number_of_ghost_cells);
    
    nested_advection.Update_Advection_Equation_Cell(grid,Z_temp,Z_forward_ghost,negative_V,dt,time+dt,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp has time n data
    
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV index=iterator.Cell_Index();
        if((*cell_mask)(index)){
            T2 Z_local=Z_forward_ghost(index)+(T).5*(Z(index)-Z_temp(index));
            if(clamp_extrema) Z_local=clamp(Z_local,Z_min(index),Z_max(index));
            Z(index)=Z_local;}
        else Z(index)=Z_forward_ghost(index);}
}
//#####################################################################
// Function Update_Advection_Equation_Face_Lookup
//#####################################################################
template<class TV,class T2,class T_NESTED_ADVECTION> void ADVECTION_MACCORMACK_UNIFORM<TV,T2,T_NESTED_ADVECTION>::
Update_Advection_Equation_Face_Lookup(const GRID<TV>& grid,ARRAY<T,FACE_INDEX<TV::m> >& Z,const ARRAY<typename TV::SCALAR,FACE_INDEX<TV::m> >& Z_ghost,const ARRAY<typename TV::SCALAR,FACE_INDEX<TV::m> >& face_velocities,
    const T dt,const T time,const ARRAY<typename TV::SCALAR,FACE_INDEX<TV::m> >* Z_min_ghost_input,const ARRAY<typename TV::SCALAR,FACE_INDEX<TV::m> >* Z_max_ghost_input,
    ARRAY<T,FACE_INDEX<TV::m> >* Z_min_input,ARRAY<T,FACE_INDEX<TV::m> >* Z_max_input,
        VECTOR<ARRAY<T2,FACE_INDEX<TV::m> >,TV::m>* flux)
{
    assert(!flux);
    assert(!ensure_second_order || !clamp_extrema);
    assert(face_mask && !Z_min_input && !Z_max_input);
    ARRAY<T,FACE_INDEX<TV::m> > Z_temp=Z;
    int number_of_ghost_cells=Z_ghost.Number_Of_Ghost_Cells();
    if(!ensure_second_order){
        ARRAY<T,FACE_INDEX<TV::m> > Z_min_ghost=Z_ghost,Z_max_ghost=Z_ghost,Z_min(grid.Cell_Indices(),no_init),Z_max(grid.Cell_Indices(),no_init);
        nested_advection.Update_Advection_Equation_Face(grid,Z_temp,Z_ghost,face_velocities,dt,time,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp now has time n+1 data
        ARRAY<T,FACE_INDEX<TV::m> > Z_forward_ghost(grid.Cell_Indices(number_of_ghost_cells),no_init);
        boundary->Fill_Ghost_Faces(grid,Z_temp,Z_forward_ghost,number_of_ghost_cells);
        boundary->Fill_Ghost_Faces(grid,Z_min,Z_min_ghost,number_of_ghost_cells);
        boundary->Fill_Ghost_Faces(grid,Z_max,Z_max_ghost,number_of_ghost_cells);
        nested_advection.Update_Advection_Equation_Face(grid,Z_temp,Z_forward_ghost,face_velocities,-dt,time+dt,&Z_min_ghost,&Z_max_ghost,&Z_min,&Z_max,0); // Z_temp has time n data
        for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
            IV index=iterator.Face_Index();int axis=iterator.Axis();
            if((*face_mask)(axis,index)){
                T Z_local=Z_forward_ghost(axis,index)+(T).5*(Z(axis,index)-Z_temp(axis,index));
                if(clamp_extrema) Z_local=clamp(Z_local,Z_min(axis,index),Z_max(axis,index));
                Z(axis,index)=Z_local;}
            else Z(axis,index)=Z_forward_ghost(axis,index);}}
    else{
        nested_advection.Update_Advection_Equation_Face(grid,Z_temp,Z_ghost,face_velocities,dt,time,0,0,0,0); // Z_temp now has time n+1 data
        ARRAY<T,FACE_INDEX<TV::m> > Z_forward_ghost(grid.Cell_Indices(number_of_ghost_cells),no_init);
        boundary->Fill_Ghost_Faces(grid,Z_temp,Z_forward_ghost,number_of_ghost_cells);
        nested_advection.Update_Advection_Equation_Face(grid,Z_temp,Z_forward_ghost,face_velocities,-dt,time+dt,0,0,0,0); // Z_temp has time n data
        for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
            FACE_INDEX<TV::m> index=iterator.Full_Index();
            if((*face_mask)(index)) Z(index)=Z_forward_ghost(index)+(T).5*(Z(index)-Z_temp(index));
            else Z(index)=Z_forward_ghost(index);}}
}

//#####################################################################
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,1>,float,ADVECTION<VECTOR<float,1>,float,ARRAY<float,FACE_INDEX<1> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,1>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<float,1>,float,AVERAGING_UNIFORM<VECTOR<float,1>,ARRAY<float,FACE_INDEX<1> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,1>,float,ARRAY<float,FACE_INDEX<1> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,1>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<float,1>,float,AVERAGING_UNIFORM<VECTOR<float,1>,ARRAY<float,FACE_INDEX<1> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,1>,float,ARRAY<float,FACE_INDEX<1> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,2>,float,ADVECTION<VECTOR<float,2>,float,ARRAY<float,FACE_INDEX<2> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,2>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<float,2>,float,AVERAGING_UNIFORM<VECTOR<float,2>,ARRAY<float,FACE_INDEX<2> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,2>,float,ARRAY<float,FACE_INDEX<2> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,2>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<float,2>,float,AVERAGING_UNIFORM<VECTOR<float,2>,ARRAY<float,FACE_INDEX<2> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,2>,float,ARRAY<float,FACE_INDEX<2> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,3>,float,ADVECTION<VECTOR<float,3>,float,ARRAY<float,FACE_INDEX<3> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,3>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<float,3>,float,AVERAGING_UNIFORM<VECTOR<float,3>,ARRAY<float,FACE_INDEX<3> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,3>,float,ARRAY<float,FACE_INDEX<3> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<float,3>,float,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<float,3>,float,AVERAGING_UNIFORM<VECTOR<float,3>,ARRAY<float,FACE_INDEX<3> > >,LINEAR_INTERPOLATION_ALL<VECTOR<float,3>,float,ARRAY<float,FACE_INDEX<3> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,1>,double,ADVECTION<VECTOR<double,1>,double,ARRAY<double,FACE_INDEX<1> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,1>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<double,1>,double,AVERAGING_UNIFORM<VECTOR<double,1>,ARRAY<double,FACE_INDEX<1> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,1>,double,ARRAY<double,FACE_INDEX<1> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,1>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<double,1>,double,AVERAGING_UNIFORM<VECTOR<double,1>,ARRAY<double,FACE_INDEX<1> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,1>,double,ARRAY<double,FACE_INDEX<1> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,2>,double,ADVECTION<VECTOR<double,2>,double,ARRAY<double,FACE_INDEX<2> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,2>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<double,2>,double,AVERAGING_UNIFORM<VECTOR<double,2>,ARRAY<double,FACE_INDEX<2> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,2>,double,ARRAY<double,FACE_INDEX<2> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,2>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<double,2>,double,AVERAGING_UNIFORM<VECTOR<double,2>,ARRAY<double,FACE_INDEX<2> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,2>,double,ARRAY<double,FACE_INDEX<2> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,3>,double,ADVECTION<VECTOR<double,3>,double,ARRAY<double,FACE_INDEX<3> > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,3>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM<VECTOR<double,3>,double,AVERAGING_UNIFORM<VECTOR<double,3>,ARRAY<double,FACE_INDEX<3> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,3>,double,ARRAY<double,FACE_INDEX<3> > > > >;
template class ADVECTION_MACCORMACK_UNIFORM<VECTOR<double,3>,double,ADVECTION_SEMI_LAGRANGIAN_UNIFORM_BETA<VECTOR<double,3>,double,AVERAGING_UNIFORM<VECTOR<double,3>,ARRAY<double,FACE_INDEX<3> > >,LINEAR_INTERPOLATION_ALL<VECTOR<double,3>,double,ARRAY<double,FACE_INDEX<3> > > > >;
}
