/*****************************************************************************
 *
 * Copyright (c) 2000 - 2006, The Regents of the University of California
 * Produced at the Lawrence Berkeley National Laboratory
 * All rights reserved.
 *
 * This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
 * full copyright notice is contained in the file COPYRIGHT located at the root
 * of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
 *
 * Redistribution  and  use  in  source  and  binary  forms,  with  or  without
 * modification, are permitted provided that the following conditions are met:
 *
 *  - Redistributions of  source code must  retain the above  copyright notice,
 *    this list of conditions and the disclaimer below.
 *  - Redistributions in binary form must reproduce the above copyright notice,
 *    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
 *    documentation and/or other materials provided with the distribution.
 *  - Neither the name of the UC/LBNL nor  the names of its contributors may be
 *    be used to endorse or promote products derived from this software without
 *    specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
 * ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
 * LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
 * DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
 * CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
 * LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
 * OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 *
 *****************************************************************************/

// ************************************************************************* //
//                            avtH5NimrodFileFormat.C                        //
// ************************************************************************* //

#include <avtH5NimrodFileFormat.h>

#include <string>

#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>

#include <Expression.h>
#include <BadIndexException.h>
#include <InvalidFilesException.h>
#include <InvalidVariableException.h>
#include <UnexpectedValueException.h>
#include <vtkStructuredGrid.h>
#include <vtkPoints.h>
#include <vtkCellType.h>

#include <vector>

#include <DebugStream.h>

using
std::string;
using
std::vector;

// ****************************************************************************
//  Function: dbg_string_attrib
//
//  Purpose:
//    Look up a string attribute from within the hdf5 file.  Echo it's value to
//    our debug logs.
//
//  Arguments:
//    id   hdf5 identifier; should be the return value from an H5GOpen call.
//    str  string to both lookup and include in the debug line.
//
//  Programmer: Tom Fogal
//  Creation:   Wed Apr 29 15:59:41 MDT 2009
//
//  Modifications:
//
// ****************************************************************************
static void dbg_string_attrib(hid_t id, const std::string &str)
{
    char *attrib = NULL;
    int h5err = H5NIMROD_read_string_attrib(id, str.c_str(), &attrib);
    if(h5err == H5NIMROD_ERR)
    {
        debug1 << "Reading '" << str << "' string in H5Nimrod gave an error!"
               << " Is this attribute guaranteed to exist for all H5Nimrod "
                  "files?" << std::endl;
    }

    // The read_string_attrib call can fail, but failure doesn't give any
    // information about whether or not the string given in the 3rd argument
    // was actually allocated.  Thus we set the string to NULL and check if it
    // got changed during the call.
    if(NULL != attrib)
    {
        debug5 << str << ": " << attrib << std::endl;
        free(attrib);
    }
}

// ****************************************************************************
//  Method: avtH5Nimrod constructor
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
//  Modifications:
//    Kathleen Bonnell, Wed Jul 2 08:43:22 PDT 2008
//    Removed unreferenced variables.
//
//    Tom Fogal, Wed Apr 29 16:11:41 MDT 2009
//    Handle errors when reading string attributes.
//
//    Jeremy Meredith, Thu Jan  7 15:36:19 EST 2010
//    Close all open ids when returning an exception.  Added error detection.
//
// ****************************************************************************

avtH5NimrodFileFormat::avtH5NimrodFileFormat (const char *filename):
    avtMTSDFileFormat (&filename, 1)
{
    // INITIALIZE DATA MEMBERS
    fname = filename;
    hid_t file;
    hid_t root_id, group_id;
    char *string_attrib;
    float time;
    file = H5Fopen (filename, H5F_ACC_RDONLY, H5P_DEFAULT);

    if (file < 0)
        EXCEPTION1 (InvalidFilesException, filename);

    hsize_t i, npoints;

    // Read attributes
    root_id = H5Gopen (file, "/");

    dbg_string_attrib(root_id, "Description");
    dbg_string_attrib(root_id, "Source");

    if (H5NIMROD_read_attrib (root_id, "time", &time) == H5NIMROD_ERR)
    {
        H5Gclose(root_id);
        H5Fclose(file);
        EXCEPTION1 (InvalidFilesException, filename);
    }

    debug5 << "time: " << time << std::endl;
    hid_t grid_id = H5Gopen (file, "/GRID");
    if (grid_id < 0)
    {
        H5Gclose(root_id);
        H5Fclose(file);
        EXCEPTION1 (InvalidFilesException, filename);
    }

    string_attrib = NULL;
    H5NIMROD_read_string_attrib (grid_id, "Coordinate System", &string_attrib);
    if(string_attrib != NULL)
    {
        debug5 << "Coordinate System: " << string_attrib << std::endl;
        if (strstr(string_attrib, "Cartesian - XYZ") == NULL)
        {
            debug5 << "Cannot handle non cartesian coordinates" << std::endl;
            H5Gclose(root_id);
            H5Gclose(grid_id);
            H5Fclose(file);
            EXCEPTION2(UnexpectedValueException, "Cartesian - XYZ",
                       string_attrib);
        }
        free(string_attrib);
    }

    string_attrib = NULL;
    H5NIMROD_read_string_attrib (grid_id, "Topology", &string_attrib);
    if (string_attrib && strstr(string_attrib, "Structured") != NULL)
    {
        structured = 1;
        debug5 << "Grid is structured!" << std::endl;
    }
    else
    {
        structured = 0;
        debug5 << "Cannot handle unstructured mesh" << std::endl;
        H5Gclose(root_id);
        H5Gclose(grid_id);
        H5Fclose (file);
        EXCEPTION2 (UnexpectedValueException, "Structured",
                    string_attrib ? string_attrib : "NULL");
    }
    if(string_attrib)
    {
        free (string_attrib);
    }

    dbg_string_attrib(grid_id, "Geometry");

    H5NIMROD_read_dims (grid_id, "X", &ndims, grid_dims);
    if (ndims != 3)
    {
        debug5 << "Cannot handle other than 3 dimensional data" << std::endl;
        H5Gclose(root_id);
        H5Gclose(grid_id);
        H5Fclose (file);
        EXCEPTION2 (UnexpectedValueException, 3, ndims);
    }
    // points
    for (i = 0, npoints = 1; i < ndims; i++)
    {
        debug5 << "ndims: " << ndims << " " << grid_dims[i] << std::endl;
        npoints *= grid_dims[i];
    }
    npoints *= ndims;
    // connectivity is implicit
    // point vars
    int num_groups = H5NIMROD_get_num_objects_matching_pattern (root_id,
            "/",
            H5G_GROUP,
            NULL);
    debug5<< "num_groups = "<<num_groups <<endl;
    char name[MAXLENGTH], name1[MAXLENGTH], name2[MAXLENGTH];

    nsteps = 0;
    for (int idx = 0; idx < num_groups; idx++)
    {
        int len_of_name = MAXLENGTH;
        H5NIMROD_get_object_name (root_id,
                "/", H5G_GROUP, idx, name, len_of_name);
        debug5<<idx<<": name= "<<name<<endl;
        if (strncmp (name, "step_", strlen ("step_")) == 0)
        {
            sprintf (name1, "/%s", name);
            stepnames.push_back( name );
            nsteps++;
            group_id = H5Gopen (root_id, name1);
            nvectorvars = H5NIMROD_get_num_objects_matching_pattern (group_id,
                    name1,
                    H5G_GROUP,
                    NULL);
            vectorvars.resize (nvectorvars);
            vectorvarnames.resize (nvectorvars);
            vectorvardims.resize (nvectorvars);
            for (int kdx = 0; kdx < nvectorvars; kdx++)
            {
                H5NIMROD_get_object_name (group_id,
                        name1,
                        H5G_GROUP, kdx, name2, len_of_name);
                vectorvarnames[kdx] = name2;
                vectorvardims[kdx] =
                    H5NIMROD_get_num_objects_matching_pattern (group_id, name2,
                            H5G_DATASET, NULL);
            }

            nscalarvars = H5NIMROD_get_num_objects_matching_pattern (group_id,
                    name1,
                    H5G_DATASET,
                    NULL);
            scalarvars.resize (nscalarvars);
            scalarvarnames.resize (nscalarvars);
            for (int kdx = 0; kdx < nscalarvars; kdx++)
            {
                H5NIMROD_get_object_name (group_id,
                        name1,
                        H5G_DATASET, kdx, name, len_of_name);
                scalarvarnames[kdx] = name;
            }
        }        
    }

    debug5 << "num scalarvars: " << nscalarvars << std::endl;
    debug5 << "num vectorvars: " << nvectorvars << std::endl;
    H5Gclose (group_id);
    H5Gclose (grid_id);
    H5Gclose (root_id);
    H5Fclose (file);
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
// ****************************************************************************

int
avtH5NimrodFileFormat::GetNTimesteps (void)
{
    return nsteps;
}


// ****************************************************************************
//  Method: avtH5NimrodFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
// ****************************************************************************

void
avtH5NimrodFileFormat::FreeUpResources (void)
{
}


// ****************************************************************************
//  Method: avtH5NimrodFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
//  Modifications:
//
//    Tom Fogal, Wed Apr 29 16:09:52 MDT 2009
//    Removed an unused variable.
//
// ****************************************************************************

void
avtH5NimrodFileFormat::PopulateDatabaseMetaData (avtDatabaseMetaData * md,
        int timeState)
{
    int nblocks = 1;                //  <-- this must be 1 for MTSD
    int block_origin = 0;
    int spatial_dimension = ndims;
    int topological_dimension = ndims;
    double *extents = NULL;


    AddMeshToMetaData (md, "Mesh", AVT_CURVILINEAR_MESH, extents, nblocks,
            block_origin, spatial_dimension, topological_dimension);

    for (int idx = 0; idx < nscalarvars; idx++)
    {
        string mesh_for_this_var = "Mesh";        // ??? -- could be multiple meshes
        string varname = scalarvarnames[idx];
        // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
        avtCentering cent = AVT_NODECENT;
        AddScalarVarToMetaData (md, varname, mesh_for_this_var, cent);
    }
    for (int idx = 0; idx < nvectorvars; idx++)
    {
        string mesh_for_this_var = "Mesh";        // ??? -- could be multiple meshes
        string varname = vectorvarnames[idx];
        // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
        avtCentering cent = AVT_NODECENT;
        AddVectorVarToMetaData (md, varname, mesh_for_this_var, cent,
                vectorvardims[idx]);
    }



}


// ****************************************************************************
//  Method: avtH5NimrodFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
// ****************************************************************************

vtkDataSet *
avtH5NimrodFileFormat::GetMesh (int timestate, const char *meshname)
{
    vtkStructuredGrid *dataset = vtkStructuredGrid::New ();
    vtkPoints *vtkpoints = vtkPoints::New ();
    vtkpoints->SetDataTypeToFloat ();
    int dims[3];

    hsize_t npoints = 1;
    for (int i = 0; i < ndims; i++)
    {
        dims[i] = (int) grid_dims[i];
        npoints *= grid_dims[i];
    }

    dataset->SetDimensions (dims);
    hid_t file;
    file = H5Fopen (fname.c_str (), H5F_ACC_RDONLY, H5P_DEFAULT);
    if (file < 0)
        EXCEPTION1 (InvalidFilesException, fname.c_str ());

    hid_t grid_id = H5Gopen (file, "/GRID");
    vtkpoints->SetNumberOfPoints (npoints);

    float *Xcoord;
    Xcoord = (float *) malloc (sizeof (float) * npoints);
    H5NIMROD_read_float32_array (grid_id, "X", NULL, ndims, NULL, Xcoord);
    _transpose_3D(Xcoord, grid_dims);
    float *Ycoord;
    Ycoord = (float *) malloc (sizeof (float) * npoints);
    H5NIMROD_read_float32_array (grid_id, "Y", NULL, ndims, NULL, Ycoord);
    _transpose_3D(Ycoord, grid_dims);
    float *Zcoord;
    Zcoord = (float *) malloc (sizeof (float) * npoints);
    H5NIMROD_read_float32_array (grid_id, "Z", NULL, ndims, NULL, Zcoord);
    _transpose_3D(Zcoord, grid_dims);

    float *pts = (float *) vtkpoints->GetVoidPointer (0);
    hsize_t idx = 0;
    for (hsize_t i = 0; i < npoints; i++)
    {
        pts[idx] = Xcoord[i];
        pts[idx + 1] = Ycoord[i];
        pts[idx + 2] = Zcoord[i];
        idx += 3;
    }
    dataset->SetPoints (vtkpoints);

    free (Xcoord);
    free (Ycoord);
    free (Zcoord);

    H5Gclose (grid_id);
    H5Fclose (file);

    return dataset;
}


// ****************************************************************************
//  Method: avtH5NimrodFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
// ****************************************************************************

vtkDataArray *
avtH5NimrodFileFormat::GetVar (int timestate, const char *varname)
{
    hid_t file;
    file = H5Fopen (fname.c_str (), H5F_ACC_RDONLY, H5P_DEFAULT);

    if (file < 0)
        EXCEPTION1 (InvalidFilesException, fname.c_str ());

    hid_t root_id, group_id;
    root_id = H5Gopen (file, "/");
    group_id = H5Gopen (root_id, stepnames[timestate].c_str ());

    float *var;
    hsize_t npoints = 1;
    for (hsize_t i = 0; i < ndims; i++)
    {
        npoints *= grid_dims[i];
    }
    var = (float *) malloc (sizeof (float) * npoints);
    H5NIMROD_read_float32_array (group_id, varname, NULL, ndims, NULL, var);
    _transpose_3D(var, grid_dims);

    hsize_t ntuples = npoints;
    vtkFloatArray *scalars = vtkFloatArray::New ();
    scalars->SetNumberOfTuples (ntuples);
    float *ptr = (float *) scalars->GetVoidPointer (0);
    memcpy (ptr, var, sizeof (float) * npoints);

    free (var);
    H5Gclose (group_id);
    H5Gclose (root_id);
    H5Fclose (file);

    return scalars;

}


// ****************************************************************************
//  Method: avtH5NimrodFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: cristina -- generated by xml2avt
//  Creation:   Fri Feb 9 08:26:27 PDT 2007
//
// ****************************************************************************

vtkDataArray *
avtH5NimrodFileFormat::GetVectorVar (int timestate, const char *varname)
{
    hid_t file;
    file = H5Fopen (fname.c_str (), H5F_ACC_RDONLY, H5P_DEFAULT);

    if (file < 0)
        EXCEPTION1 (InvalidFilesException, fname.c_str ());

    hid_t root_id, group_id;
    root_id = H5Gopen (file, "/");
    group_id = H5Gopen (root_id, stepnames[timestate].c_str ());
    int num_comp = H5NIMROD_get_num_objects_matching_pattern (group_id,
            varname,
            H5G_DATASET,
            NULL);

    hsize_t npoints = 1;
    for (hsize_t i = 0; i < ndims; i++)
    {
        npoints *= grid_dims[i];
    }

    float **comp;
    comp = (float **) malloc (sizeof (float *) * num_comp);
    for (int idx = 0; idx < num_comp; idx++)
        comp[idx] = (float *) malloc (sizeof (float) * npoints);

    hid_t vector_id = H5Gopen (group_id, varname);
    char name[MAXLENGTH];
    int len_of_name = MAXLENGTH;
    for (int idx = 0; idx < num_comp; idx++)
    {
        H5NIMROD_get_object_name (group_id,
                varname, H5G_DATASET, idx, name, len_of_name);
        H5NIMROD_read_float32_array (vector_id,
                name, NULL, ndims, NULL, comp[idx]);
        _transpose_3D(comp[idx], grid_dims);
    }

    vtkFloatArray *vector = vtkFloatArray::New ();
    hsize_t ntuples = npoints;
    vector->SetNumberOfComponents (num_comp);
    vector->SetNumberOfTuples (ntuples);
    float *ptr = (float *) vector->GetVoidPointer (0);
    hsize_t count = 0;
    for (hsize_t idx = 0; idx < ntuples; idx++)
    {
        for (int jdx = 0; jdx < num_comp; jdx++)
        {
            ptr[count + jdx] = comp[jdx][idx];
        }
        count += num_comp;
    }
    for (int idx = 0; idx < num_comp; idx++)
    {
        if (comp[idx] != NULL)
            free (comp[idx]);
    }
    if (comp != NULL)
        free (comp);

    return vector;

    H5Gclose (vector_id);
    H5Gclose (group_id);
    H5Gclose (root_id);
    H5Fclose (file);
}