#include <Core/Grids/GRID.h>
#include <Core/Parsing/PARSE_ARGS.h>
#include <Core/Particles/PARTICLES.h>
#include <Core/Read_Write/FILE_UTILITIES.h>
#include <fstream>
#include <iostream>
#include <regex>

using namespace PhysBAM;

// usage ./dump_grid  output/x/node_velocities [-3d]

template<class T,int dim> void read_node_velocity(std::string input_directory)
{
    typedef VECTOR<T,dim> TV;
    typedef VECTOR<int,dim> IV;

    ARRAY<TV,IV> V;
    Read_From_File(input_directory,V);

    GRID<TV> grid;
    GRID<TV> mac_grid;
    std::string common_grid_file,node_velocities_file;
    // read directory and frame number
    std::regex r("[0-9]+/node_velocities(\\.gz)?");
    common_grid_file=regex_replace(input_directory,r,"common/grid");
    LOG::cout<<input_directory<<" "<<common_grid_file<<std::endl;
    if(File_Exists(common_grid_file)){
        LOG::cout<<"Reading grid from '"<<common_grid_file<<"'..."<<std::endl;
        Read_From_File(common_grid_file,mac_grid);
        grid=mac_grid.Get_Regular_Grid_At_MAC_Positions();
    }

    ARRAY<TV> vector_field,vector_locations;
    int idx=0;
    vector_field.Resize(V.Size().Product());
    vector_locations.Resize(V.Size().Product());
    // TODO only get 2d here; don't have 3d see -> OPENGL_GRID_BASED_VECTOR_FIELD_3D->Update()
    for(RANGE_ITERATOR<TV::m> it(V.domain);it.Valid();it.Next()){
      vector_field(idx)=V(it.index);
      vector_locations(idx)=grid.Node_X(it.index);
      idx++;}
    // vector_field + vector_locations will do the drawing
    // LOG::printf("mac_grid dX %P grid dX %P \n",mac_grid.dX,grid.dX);
    LOG::printf("vector field %P \nvector location %P\n",vector_field,vector_locations);
    // LOG::printf("%P max xize %P\n",V.Size(),V.Pack_Size());
    // LOG::printf("V is %P\n",V);
}

int main(int argc,char *argv[])
{
    typedef double T;
    std::string input_directory;
    bool use_3d=false;
    PARSE_ARGS parse_args(argc,argv);
    parse_args.Extra(&input_directory,"file","grid data");
    parse_args.Add("-3d",&use_3d,"run 3D examples");
    parse_args.Parse();

    if(!use_3d)
        read_node_velocity<T,2>(input_directory);
    else
        read_node_velocity<T,3>(input_directory);

    return 0;
}

// $PHYSBAM/Tools/dump_particle/dump_grid output/2/node_velocities.gz 
