//#####################################################################
// Copyright 2009-2010, Michael Lentine, 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_DEFINITION.h>
#include <Core/Geometry_Particles/DEBUG_PARTICLES.h>
#include <Core/Interpolation/AVERAGING_UNIFORM.h>
#include <Core/Interpolation/FACE_LOOKUP_UNIFORM.h>
#include <Core/Interpolation/LINEAR_INTERPOLATION.h>
#include <Core/Log/LOG.h>
#include <Core/Read_Write/FILE_UTILITIES.h>
#include "SMOKE_EXAMPLE.h"
using namespace PhysBAM;
//#####################################################################
// Constructor
//#####################################################################
template<class TV> SMOKE_EXAMPLE<TV>::
SMOKE_EXAMPLE(const STREAM_TYPE stream_type_input,int number_of_threads)
    :stream_type(stream_type_input),
    debug_particles(DEBUG_PARTICLES<TV>::Get()),ghost(5),
    last_frame(100),frame_dt((T)1/24),
    restart(0),write_debug_data(true),viewer_dir("output"),N_boundary(false),
    debug_divergence(false),alpha(0.1),beta(0.00366),
    cfl(.9),grid(IV(),RANGE<TV>::Unit_Box()),
    projection(grid,false,false),boundary(0)
{
    for(int i=0;i<TV::m;i++){domain_boundary(i)(0)=true;domain_boundary(i)(1)=true;}
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> SMOKE_EXAMPLE<TV>::
~SMOKE_EXAMPLE()
{
}
//#####################################################################
// Function CFL
//#####################################################################
template<class TV> typename TV::SCALAR SMOKE_EXAMPLE<TV>::
CFL(ARRAY<T,FACE_INDEX<TV::m> >& face_velocities)
{
    T dt=FLT_MAX;
    T dt_convection=0;
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        IV cell=iterator.Cell_Index();T local_V_norm=0;
        for(int axis=0;axis<GRID<TV>::dimension;axis++)
            local_V_norm+=grid.one_over_dX[axis]*maxabs(face_velocities(axis,grid.First_Face_Index_In_Cell(axis,cell)),face_velocities(axis,grid.Second_Face_Index_In_Cell(axis,cell)));
        dt_convection=max(dt_convection,local_V_norm);}
    dt=min(dt,(T)1.0/dt_convection);
    return dt;
}
//#####################################################################
// Function Time_At_Frame
//#####################################################################
template<class TV> typename TV::SCALAR SMOKE_EXAMPLE<TV>::
Time_At_Frame(const int frame) const 
{
    return frame*frame_dt;
}
//#####################################################################
// Function Initialize_Grid
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Initialize_Grid(IV counts,RANGE<TV> domain) 
{
    grid.Initialize(counts,domain);
}
//#####################################################################
// Function Initialize_Fields
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Initialize_Fields() 
{
    for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()) face_velocities(iterator.Full_Index())=0;
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()) density(iterator.Cell_Index())=0;
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()) temperature(iterator.Cell_Index())=273; // add temperature
}
//#####################################################################
// Function Get_Scalar_Field_Sources
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Get_Scalar_Field_Sources(const T time)
{
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next())
    {
        if(source1.Lazy_Inside(iterator.Location())) {density(iterator.Cell_Index())=1;}
        if(source1.Lazy_Inside(iterator.Location())) {temperature(iterator.Cell_Index())=275;}    
        //else if(source2.Lazy_Inside(iterator.Location())) {density(iterator.Cell_Index())=1;}
    }
}
//#####################################################################
// Function Set_Boundary_Conditions
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Set_Boundary_Conditions(const T time, const T source_velocities)
{
    projection.elliptic_solver->psi_D.Fill(false);projection.elliptic_solver->psi_N.Fill(false);
    for(int axis=0;axis<TV::m;axis++) for(int axis_side=0;axis_side<2;axis_side++){int side=2*axis+axis_side;
        if(domain_boundary(axis)(axis_side)){
            IV interior_cell_offset=axis_side==0?IV():-IV::Axis_Vector(axis);    
            for(FACE_ITERATOR<TV> iterator(grid,1,BOUNDARY_REGION,side);iterator.Valid();iterator.Next()){IV cell=iterator.Face_Index()+interior_cell_offset;
                IV boundary_face=axis_side==0?iterator.Face_Index()+IV::Axis_Vector(axis):iterator.Face_Index()-IV::Axis_Vector(axis);
                if(N_boundary){
                    FACE_INDEX<TV::m> face(axis,boundary_face);
                    projection.elliptic_solver->psi_N(face)=true;
                    face_velocities(face)=0;}
                else projection.elliptic_solver->psi_D(cell)=true;
                projection.p(cell)=0;}}}
    for(FACE_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){
        if(source1.Lazy_Inside(iterator.Location())){
            projection.elliptic_solver->psi_N(iterator.Full_Index())=true;
            if(iterator.Axis()==1)face_velocities(iterator.Full_Index())=source_velocities;
            else face_velocities(iterator.Full_Index())=0;}
        // else if(source2.Lazy_Inside(iterator.Location())){   // second source
        //     projection.elliptic_solver->psi_N(iterator.Full_Index())=true;
        //     if(iterator.Axis()==1)face_velocities(iterator.Full_Index())=-source_velocities;
        //     else face_velocities(iterator.Full_Index())=0;}
    }
}
//#####################################################################
// Function Register_Output_Files
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Register_Output_Files()
{
    if(TV::m==3) viewer_output.template Add_Axes<TV::m>("axes",{});
    viewer_output.Add("grid",grid,true,{.draw=(TV::m==2)});
    viewer_output.Add("mac_velocities","grid",face_velocities,false,{.draw=false});
    viewer_output.Add("debug_particles",DEBUG_PARTICLES<TV>::Get(),false,{.draw=false});
    viewer_output.Add("density","grid",density,false,{.draw=false});
    viewer_output.Add("temperature","grid",temperature,false,{.draw=false});
    if(write_debug_data)
    {
        viewer_output.Add("pressure","grid",projection.p,false,{.draw=false});
        viewer_output.Add("psi_N","grid",projection.elliptic_solver->psi_N,false,{.draw=false});
        viewer_output.Add("psi_D","grid",projection.elliptic_solver->psi_D,false,{.draw=false});
    }
}
//#####################################################################
// Function Read_Output_Files
//#####################################################################
template<class TV> void SMOKE_EXAMPLE<TV>::
Read_Output_Files()
{
    Read_From_File(viewer_dir.Frame_File("density"),density);
    std::string filename;
    filename=viewer_dir.Frame_File("mac_velocities");
    if(File_Exists(filename)){
        LOG::cout<<"Reading mac_velocities "<<filename<<std::endl;
        Read_From_File(filename,face_velocities);}
    filename=viewer_dir.Frame_File("pressure");
    if(File_Exists(filename)){
        LOG::cout<<"Reading pressure "<<filename<<std::endl;
        Read_From_File(filename,projection.p);}
}
//#####################################################################
namespace PhysBAM{
template class SMOKE_EXAMPLE<VECTOR<float,2> >;
template class SMOKE_EXAMPLE<VECTOR<float,3> >;
template class SMOKE_EXAMPLE<VECTOR<double,2> >;
template class SMOKE_EXAMPLE<VECTOR<double,3> >;
}
