/*****************************************************************************
*
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* 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 LLNS/LLNL nor the names of  its contributors may
*    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.
*
*****************************************************************************/

#include <vector>
#include <snprintf.h>
#include <netcdf.h>

#include <avtBasicNETCDFReader.h>
#include <NETCDFFileObject.h>
#include <avtDatabase.h>
#include <avtDatabaseMetaData.h>
#include <avtSTSDFileFormatInterface.h>
#include <DebugStream.h>

#include <vtkCellArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkPoints.h>

#include <vtkUnsignedCharArray.h>
#include <vtkShortArray.h>
#include <vtkIntArray.h>
#include <vtkLongArray.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkVisItUtility.h>

#include <InvalidVariableException.h>

#define TIME_DIMENSION      -1

#ifdef PARALLEL
#include <avtParallel.h>
#endif

// ****************************************************************************
// Method: avtBasicNETCDFReader::avtBasicNETCDFReader
//
// Purpose: 
//   Constructor for the avtBasicNETCDFReader class.
//
// Arguments:
//   filename : The name of the file being read.
//   f        : The file object associated with the file being read.
//
// Programmer: Brad Whitlock
// Creation:   Thu Aug 18 18:03:59 PST 2005
//
// Modifications:
//   Brad Whitlock, Wed Apr 26 17:40:20 PST 2006
//   Initialized meshNamesCreated.
//
//   Mark C. Miller, Tue Aug 15 15:28:11 PDT 2006
//   Added procNum, procCount to support on-the-fly parallel decomposition
//
// ****************************************************************************

avtBasicNETCDFReader::avtBasicNETCDFReader(const char *filename) :
    avtNETCDFReaderBase(filename), meshNameToDimensionsSizes(), meshNameToNCDimensions(), varToDimensionsSizes()
{
    meshNamesCreated = false;
    procNum = 0;
    procCount = 1;
#ifdef PARALLEL
    procNum = PAR_Rank();
    procCount = PAR_Size();
#endif
}

avtBasicNETCDFReader::avtBasicNETCDFReader(const char *filename, NETCDFFileObject *f) :
    avtNETCDFReaderBase(filename, f), meshNameToDimensionsSizes(), meshNameToNCDimensions(), varToDimensionsSizes()
{
    meshNamesCreated = false;
    procNum = 0;
    procCount = 1;
#ifdef PARALLEL
    procNum = PAR_Rank();
    procCount = PAR_Size();
#endif
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::~avtBasicNETCDFReader
//
// Purpose: 
//   Destructor for the avtBasicNETCDFReader class.
//
// Programmer: Brad Whitlock
// Creation:   Thu Aug 18 18:04:36 PST 2005
//
// Modifications:
//   
// ****************************************************************************

avtBasicNETCDFReader::~avtBasicNETCDFReader()
{
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::CreateGlobalAttributesString
//
// Purpose: 
//   Create a string for the global attributes.
//
// Arguments:
//   nGlobalAtts : The number of global attributes
//   gaString    : global attributes string.
//
// Returns:    
//
// Note:       
//
// Programmer: Brad Whitlock
// Creation:   Fri Oct 30 11:01:08 PDT 2009
//
// Modifications:
//   
// ****************************************************************************

void
avtBasicNETCDFReader::CreateGlobalAttributesString(int nGlobalAtts, std::string &gaString)
{
    for(int i = 0; i < nGlobalAtts; ++i)
    {
        int     status;
        char    attname[NC_MAX_NAME+1];
        nc_type atttype;
        size_t  attsize;
        if((status = nc_inq_attname(fileObject->GetFileHandle(), NC_GLOBAL, i, attname))
            == NC_NOERR)
        {
            if((status = nc_inq_att(fileObject->GetFileHandle(), NC_GLOBAL, attname, &atttype,
                                    &attsize)) == NC_NOERR)
            {
                std::string tmpStr("\t");
                tmpStr += attname;

                if(atttype == NC_CHAR)
                {
                    char *value = new char[attsize+1];
                    nc_get_att_text(fileObject->GetFileHandle(), NC_GLOBAL, attname, value);
                    value[attsize] = '\0';
                    char *c2 = value + attsize - 1;
                    while(c2 >= value && *c2 == ' ')
                        *c2-- = '\0';
                    tmpStr += " = \"";
                    tmpStr += value;
                    tmpStr += "\"";

                    delete [] value;
                }

                tmpStr += "\n";
                gaString += tmpStr;
            }
        }
    }
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::PopulateDatabaseMetaData
//
// Purpose: 
//   Populates the metadata from information in the file.
//
// Arguments:
//   md : The metadata object to populate.
//
// Programmer: Brad Whitlock
// Creation:   Thu Aug 18 18:05:29 PST 2005
//
// Modifications:
//    Jeremy Meredith, Thu Aug 25 12:55:29 PDT 2005
//    Added group origin to mesh metadata constructor.
//
//    Brad Whitlock, Wed Apr 26 17:53:24 PST 2006
//    I made it possible to call with a NULL metadata pointer so we can call
//    this method on ActivateTimestep to ensure that the 
//    meshNameToDimensionsSizes map is populated. I also added a check to set 
//    the validVariable flag for meshes and variables so that we don't have 
//    problems later. By adding the valid variable check, I was able to make 
//    4+ dimensional arrays safely appear in the variable lists.
//
//    Mark C. Miller, Tue Aug 15 15:28:11 PDT 2006
//    Added call to SetFormatCanDoDomainDecomposition
//
//    Jeremy Meredith, Mon Aug 28 17:39:33 EDT 2006
//    Added test to make sure a variable had >0 dimensions before examining
//    them.
//
//    Hank Childs, Fri Feb 23 09:20:38 PST 2007
//    Fix memory leak.
//
//    Brad Whitlock, Thu Jul 15 12:15:56 PDT 2010
//    I added support for curves over time.
//
//    Kathleen Bonnell, Wed Sep 8 20:31:55 MST 2010 
//    Check value of maxDimNotOne before using it as an index, to prevent SEGV.
//
//    Brad Whitlock, Wed Jan  4 15:47:43 PST 2012
//    Handle missing data conventions.
//
//    Mark C. Miller, Wed Nov 20 14:42:56 PST 2013
//    Two passes to add original (node-centered) variables and also an
//    "as_zonal/" submenu for the same variables handled as zone centered. 
// ****************************************************************************

void
avtBasicNETCDFReader::PopulateDatabaseMetaData(int timeState, avtDatabaseMetaData *md)
{
    debug4 << "avtBasicNETCDFReader::PopulateDatabaseMetaData" << endl;
    if(DebugStream::Level4())
        fileObject->PrintFileContents(DebugStream::Stream4());

    int status, nDims, nVars, nGlobalAtts, unlimitedDimension;
    status = nc_inq(fileObject->GetFileHandle(), &nDims, &nVars, &nGlobalAtts,
                    &unlimitedDimension);
    debug4 << "nDims=" << nDims
           << ", nVars=" << nVars
           << ", nGlobalAtts=" << nGlobalAtts
           << ", unlimitedDimension=" << unlimitedDimension
           << endl;
    if(status != NC_NOERR)
    {
        fileObject->HandleError(status);
        return;
    }

    if(md != 0)
    {
        std::string gaString;
        CreateGlobalAttributesString(nGlobalAtts, gaString);
        md->SetDatabaseComment(gaString);
    }

    // Get the size of all of the dimensions in the file.
    size_t *dimSizes = new size_t[nDims];
    int i;
    for(i = 0; i < nDims; ++i)
    {
        char   dimName[NC_MAX_NAME+1];
        size_t dimSize;
        if((status = nc_inq_dim(fileObject->GetFileHandle(), i, dimName,
           &dimSize)) == NC_NOERR)
        {
            dimSizes[i] = dimSize;
        }
        else
        {
            dimSizes[i] = 1;
            fileObject->HandleError(status);
        }
    }

    // Get the time dimension.
    int timedim = -1, time_nts = 0;
    std::string timedimname;
    GetTimeDimension(fileObject, timedim, time_nts, timedimname);

    meshNameToDimensionsSizes.clear();
    meshNameToNCDimensions.clear();
    varToDimensionsSizes.clear();

    // Iterate over the variables and create a list of meshes names and add
    // the variable to the metadata.
    for(i = 0; i < nVars; ++i)
    {
        char varname[NC_MAX_NAME+1];
        nc_type vartype;
        int  varndims;
        int  vardims[NC_MAX_VAR_DIMS];
        int  varnatts;
        if((status = nc_inq_var(fileObject->GetFileHandle(), i, varname,
                                &vartype, &varndims, 
                                vardims, &varnatts)) == NC_NOERR)
        {
            // Count the number of variable dimensions that are greater
            // than 1.
            if (varndims == 0)
                continue;

            int nGt1Dims = 0;
            int maxDim = dimSizes[vardims[0]];
            int maxDimIndex = 0;
            int maxNonTimeDim = 0;
            int maxNonTimeDimIndex = -1;
            bool hasTimeDimension = false;
            for(int dim = 0; dim < varndims; ++dim)
            {
                int d = dimSizes[vardims[dim]];
                if(d > 1)
                    ++nGt1Dims;
                if(d > maxDim)
                {
                    maxDim = d;
                    maxDimIndex = dim;
                }

                if(vardims[dim] == timedim)
                    hasTimeDimension = true;
                else if(d > 1 && (maxNonTimeDimIndex == -1 || d > maxNonTimeDim))
                {
                    maxNonTimeDim = d;
                    maxNonTimeDimIndex = dim;
                }
            }

            if(nGt1Dims == 1 || (varndims == 2 && nGt1Dims == 2 && hasTimeDimension))
            {
                if(maxDim > 1 && 
                   (vartype == NC_INT ||
                    vartype == NC_SHORT||
                    vartype == NC_FLOAT ||
                    vartype == NC_DOUBLE)
                  )
                {
                    int xdim = maxDimIndex;
                    int xdimSize = maxDim;

                    // Time varying curve.
                    if(nGt1Dims == 2 && hasTimeDimension)
                    {
                        xdim = maxNonTimeDimIndex;
                        xdimSize = maxNonTimeDim;

                        // We only save the xdim for time varying curves so we can
                        // tell them apart from normal curves.
                        intVector meshDims; meshDims.push_back(xdim);
                        meshNameToNCDimensions[varname] = meshDims;

                        intVector vDims;
                        for(int dim = varndims-1; dim >= 0; --dim)
                        {
                            int d;
                            if(vardims[dim] != timedim)
                                d = (int)dimSizes[vardims[dim]];
                            else
                                d = TIME_DIMENSION;
                            vDims.push_back(d);
                        }
                        varToDimensionsSizes[varname] = vDims;
                    }

                    if(md != 0)
                    {
                        avtCurveMetaData *cmd = new avtCurveMetaData;
                        cmd->name = varname;
                        cmd->yLabel = varname;
                        fileObject->ReadStringAttribute(varname, "units",
                            cmd->yUnits);
                       
                        char dimName[NC_MAX_NAME+1];
                        status = nc_inq_dimname(fileObject->GetFileHandle(),
                                                vardims[xdim], dimName);
                        if(status == NC_NOERR)
                        {
                            cmd->xLabel = dimName;
                            fileObject->ReadStringAttribute(dimName, "units",
                                                            cmd->xUnits);
                        }
                        md->Add(cmd);
                    }

                    intVector meshDimsSizes; meshDimsSizes.push_back(xdimSize);
                    meshNameToDimensionsSizes[varname] = meshDimsSizes;
                }
            }
            else
            {
                char tmp[100];
                std::string meshName("mesh");
                intVector meshDims, meshDimSizes, vDims;
                int elems = 1;
                for(int dim = varndims-1; dim >= 0; --dim)
                {
                    int d;
                    if(vardims[dim] != timedim)
                    {
                        vDims.push_back(vardims[dim]);
                        d = (int)dimSizes[vardims[dim]];
                        elems *= d;
                    }
                    else
                        d = TIME_DIMENSION;

                    meshDims.push_back(vardims[dim]);
                    meshDimSizes.push_back(d);
                }

                int maxDimNotOne = meshDimSizes.size();
                while(maxDimNotOne >= 1 && meshDimSizes[maxDimNotOne-1] == 1)
                     maxDimNotOne--;

                // Make sure it's not 1 or 0 in some of its dimensions
                if(meshDimSizes.size() >= 2 && elems != maxDim && elems > 0)
                {
                    // Come up with the mesh name that we'll use for this
                    // variable. Note that we don't consider the time dimension
                    // so time-varying and non-time-varying variables can share
                    // the same dimension if possible.
                    int nSpatialDims = 0;
                    for(int j = 0; j < maxDimNotOne; ++j)
                    {
                        if(meshDimSizes[j] != TIME_DIMENSION)
                        {
                            if(j != 0)
                                SNPRINTF(tmp, 100, "x%d", meshDimSizes[j]);
                            else
                                SNPRINTF(tmp, 100, "%d", meshDimSizes[j]);
                            ++nSpatialDims;
                            meshName += tmp;
                        }
                    }

                    // Add the name of the mesh to the list of meshes.
                    if(meshNameToDimensionsSizes.find(meshName) == meshNameToDimensionsSizes.end())
                    { 
                        // Filter out time from the dimensions so time varying and static
                        // variables can share the same mesh.
                        intVector meshDimsWithoutTime, meshDimSizesWithoutTime;
                        for(size_t dim = 0; dim < meshDimSizes.size(); ++dim)
                            if(meshDimSizes[dim] != TIME_DIMENSION)
                            {
                                meshDimSizesWithoutTime.push_back(meshDimSizes[dim]);
                                meshDimsWithoutTime.push_back(meshDims[dim]);
                            }
                        meshNameToDimensionsSizes[meshName] = meshDimSizesWithoutTime;
                        meshNameToNCDimensions[meshName] = meshDimsWithoutTime;

                        if(md != 0)
                        {
                            md->SetFormatCanDoDomainDecomposition(true);
                            for (int pass = 0; pass < 2; pass++)
                            {
                                avtMeshMetaData *mmd = new avtMeshMetaData(
                                    (pass?meshName:("as_zonal/"+meshName)), 
                                    1, 1, 1, 0, nSpatialDims, nSpatialDims,
                                    AVT_RECTILINEAR_MESH);
                                mmd->validVariable = nSpatialDims <= 3;
                                // Provide the dimension names as the axis labels.
                                char dimName[NC_MAX_NAME+1];
                                size_t sz;
                                if((status = nc_inq_dim(fileObject->GetFileHandle(), vDims[0], dimName,
                                    &sz)) == NC_NOERR)
                                {
                                    mmd->xLabel = dimName;
                                    fileObject->ReadStringAttribute(dimName, "units", mmd->xUnits);
                                }
                                if((status = nc_inq_dim(fileObject->GetFileHandle(), vDims[1], dimName,
                                    &sz)) == NC_NOERR)
                                {
                                    mmd->yLabel = dimName;
                                    fileObject->ReadStringAttribute(dimName, "units", mmd->yUnits);
                                }
                                if(nSpatialDims == 3)
                                {
                                    if((status = nc_inq_dim(fileObject->GetFileHandle(), vDims[2], dimName,
                                        &sz)) == NC_NOERR)
                                    {
                                        mmd->zLabel = dimName;
                                        fileObject->ReadStringAttribute(dimName, "units", mmd->zUnits);
                                    }
                                }
                                md->Add(mmd);
                            }
                        }
                    }

                    // Try and get the variable units.
                    varToDimensionsSizes[varname] = meshDimSizes;
                    if(md != 0)
                    {
                        for (int pass = 0; pass < 2; pass++)
                        {
                            avtScalarMetaData *smd = new avtScalarMetaData(
                                (pass?varname:("as_zonal/"+std::string(varname))),
                                (pass?meshName:("as_zonal/"+meshName)),
                                pass?AVT_NODECENT:AVT_ZONECENT);
                            smd->hasUnits = fileObject->ReadStringAttribute(
                                varname, "units", smd->units);
                            smd->validVariable = nSpatialDims <= 3;
                            HandleMissingData(varname, smd);
                            md->Add(smd);
                        }
                    }
                } 
            }
        }
    }

    delete [] dimSizes;
    meshNamesCreated = true;
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::ReturnSpatialDimensionIndices
//
// Purpose: 
//   Returns the indices of the spatial dimensions.
//
// Programmer: Brad Whitlock
// Creation:   Thu Oct 29 15:41:42 PDT 2009
//
// Modifications:
//   
// ****************************************************************************

bool
avtBasicNETCDFReader::ReturnSpatialDimensionIndices(const intVector &dims, int sDims[3], int &nSDims) const
{
    size_t i;
    const char *mName = "avtBasicNETCDFReader::ReturnValidDimensions: ";

    // Look for up to 3 valid spatial dimensions.
    nSDims = 0;
    for(i = 0; i < dims.size() && nSDims < 3; ++i)
    {
        if(dims[i] > 1 && dims[i] != TIME_DIMENSION)
        {
            sDims[nSDims++] = i;
        }
    }

    // Count the number of cells that comprise the spatial dimensions
    int nCells = 1;
    debug5 << mName << "validDims=(";
    for(i = 0; i < (size_t)nSDims; ++i)
    {
        nCells *= dims[sDims[i]];
        debug5 << dims[sDims[i]] << ", ";
    }
    debug5 << ")" << endl;

    // Count the number of cells that comprise all dimensions (except time).
    int nValues = 1;
    debug5 << mName << "actualDims=(";
    for(i = 0; i < dims.size(); ++i)
    {
        if(dims[i] != TIME_DIMENSION)
            nValues *= dims[i];
        debug5 << dims[i] << ", ";
    }
    debug5 << ")" << endl;

    return nCells == nValues;
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::ReturnDimStartsAndCounts
//
// Purpose: 
//   Returns the dimStarts and dimCounts for all dimensions.
//
// Programmer: Brad Whitlock
// Creation:   Thu Oct 29 15:41:42 PDT 2009
//
// Modifications:
//   Mark C. Miller, Tue Aug 15 15:28:11 PDT 2006
//   Added code to support on-the-fly domain decomposition 
//
//   Mark C. Miller, Wed Aug 16 14:45:22 PDT 2006
//   Fixed possible exit without initializing all contents of starts/counts
//
//   Mark C. Miller, Tue Dec  5 18:14:58 PST 2006
//   Fixed UMRs
//
// ****************************************************************************

void
avtBasicNETCDFReader::ReturnDimStartsAndCounts(int timeState, const intVector &dims, 
    intVector &dimStarts, intVector &dimCounts) const
{
    //
    // Initialize starts/counts, slicing in time.
    //
    dimStarts.clear();
    dimCounts.clear();
    for (size_t i = 0; i < dims.size(); i++)
    {
        if(dims[i] == TIME_DIMENSION)
        {
            dimStarts.push_back(timeState);
            dimCounts.push_back(1);
        }
        else
        {
            dimStarts.push_back(0);
            dimCounts.push_back(dims[i]);
        }
    }

    // Get the indices of the spatial dimensions
    int spatialDimIndices[3] = {0,0,0}, nSpatialDims = 0;
    ReturnSpatialDimensionIndices(dims, spatialDimIndices, nSpatialDims);

    // Compute how many values that makes.
    int spatialSizes[3] = {0,0,0};
    int nValues = 1;
    for(int i = 0; i < nSpatialDims; ++i)
    {
        spatialSizes[i] = dims[spatialDimIndices[i]];
        nValues *= spatialSizes[i];
    }

    //
    // We won't decompose something that is smaller than some threshold
    //
    if (nValues < 100000)
        return;

    //
    // Here's where we let VisIt decide how to divide up the spatial dimensions.
    //

    //
    // Above, we're counting nodes (e.g. mesh lines). The decomposition
    // stuff operates on zones.
    //
    int validZDims[3] = {0,0,0};
    for (int i = 0; i < nSpatialDims; i++)
        validZDims[i] = spatialSizes[i]-1;

    //
    // Ok, now compute the zone-oriented domain decomposition
    //
    int domCount[3] = {0, 0, 0};
    avtDatabase::ComputeRectilinearDecomposition(
        nSpatialDims, procCount,
        validZDims[0], validZDims[1], validZDims[2],
        &domCount[0], &domCount[1], &domCount[2]);

    debug5 << "Decomposition: " << domCount[0] << ", "
           << domCount[1] << ", " << domCount[2] << endl;

    //
    // Determine this processor's logical domain (e.g. domain ijk) indices
    //
    int domLogicalCoords[3] = {0, 0, 0};
    avtDatabase::ComputeDomainLogicalCoords(nSpatialDims, domCount, procNum,
        domLogicalCoords);

    debug5 << "Processor " << procNum << " domain logical coords: "
           << domLogicalCoords[0] << ", " << domLogicalCoords[1] << ", "
           << domLogicalCoords[2] << endl;

    //
    // compute the bounds, in terms of output zone numbers,
    // of this processor's domain. Store the divided spatial domain's starts
    // and counts, overwriting the previous values.
    //
    debug5 << "Processor " << procNum << " zone-centered bounds..." << endl;
    for (int i = 0; i < nSpatialDims; i++)
    {
        avtDatabase::ComputeDomainBounds(validZDims[i], domCount[i], domLogicalCoords[i],
            &dimStarts[spatialDimIndices[i]], &dimCounts[spatialDimIndices[i]]);
        dimCounts[spatialDimIndices[i]]++; // convert to # of zones to # of nodes  
        debug5 << "   start[" << i << "] = " << dimStarts[spatialDimIndices[i]]
               << ",  count[" << i << "] = " << dimCounts[spatialDimIndices[i]]-1 << endl;
    }
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::GetMesh
//
// Purpose: 
//   Returns the specified mesh.
//
// Arguments:
//   timeState : The time state.
//   var       : The name of the mesh to create.
//
// Returns:    A vtkDataSet containing the mesh or 0.
// 
// Programmer: Brad Whitlock
// Creation:   Thu Aug 18 18:05:59 PST 2005
//
// Modifications:
//   Brad Whitlock, Wed Apr 26 17:56:15 PST 2006
//   I made it call PopulateDatabaseMetaData.
//
//   Mark C. Miller, Tue Aug 15 15:28:11 PDT 2006
//   Added code to support on-the-fly domain decomp
//
//   Kathleen Bonnell, Mon Jul 14 14:07:39 PDT 2008
//   Specify curves as 1D rectilinear grids with yvalues stored in point data.
//
//   Brad Whitlock, Thu Oct 29 15:38:29 PDT 2009
//   I rewrote it with support for time.
//
//   Brad Whitlock, Thu Jul 15 12:13:30 PDT 2010
//   I added support for curves over time.
//
//   Mark C. Miller, Wed Nov 20 14:42:56 PST 2013
//   Handle mesh construction for zone centered variants of variables too. 
//
//   Burlen Loring, Mon Jul 14 13:03:55 PDT 2014
//   fix out of bounds index into vector of dimension id's. added a check
//   to prevent out of bounds access into coordvals array (x-axis values)
//   when it's a different length than the variable (y-axis values). In
//   that case id values are used instead.
//
//   Kathleen Biagas, Thu Jun 4 18:18:42 MST 2015
//   Create base_index FieldData array.
//
// ****************************************************************************

vtkDataSet *
avtBasicNETCDFReader::GetMesh(int timeState, const char *var)
{
    const char *mName = "avtBasicNETCDFReader::GetMesh: ";
    debug4 << mName << "var=" << var << endl;
    vtkDataSet *retval = 0;
    char const *bare_var = var;
    if (!strncmp(var, "as_zonal/", sizeof("as_zonal/")-1))
        bare_var = &var[sizeof("as_zonal/")-1];

    avtCentering centering = AVT_NODECENT;
    if (!strncmp(var, "as_zonal/", sizeof("as_zonal/")-1))
        centering = AVT_ZONECENT;
 
    // Populate the mesh names if we've not done so yet.
    if(!meshNamesCreated)
        PopulateDatabaseMetaData(timeState, 0);

    StringIntVectorMap::const_iterator mesh = meshNameToDimensionsSizes.find(bare_var);
    if(mesh != meshNameToDimensionsSizes.end())
    {
        if(mesh->second.size() == 1)
        {
            StringIntVectorMap::const_iterator meshDims = meshNameToNCDimensions.find(bare_var);
            if(meshDims != meshNameToNCDimensions.end())
            {
                // Time varying curves

                // Try and get values for the x-axis from a dimension. If we can't 
                // get it then we just use indices.
                float *coordvals = 0;                
                int status;
                size_t sz = 0;
                char dimName[NC_MAX_NAME+1];
                debug4 << mName << "Looking for X coordinate dimension name" << endl;
                if((status = nc_inq_dim(fileObject->GetFileHandle(),
                    meshDims->second[0], dimName, &sz)) == NC_NOERR)
                {
                    TypeEnum dvart = NO_TYPE;
                    int dvarndims = 0;
                    int *dvardims = 0;
                    debug4 << mName << "Looking for " << dimName << " array as X coordinate" << endl;
                    if(fileObject->InqVariable(dimName, &dvart, &dvarndims, &dvardims))
                    {
                        if(dvarndims == 1 && (size_t)dvardims[0] == sz)
                            coordvals = ReadArray(dimName);
                        delete [] dvardims;
                    }
                }
 
                // Read the y value using GetVar
                vtkDataArray *yv = GetVar(timeState, var);
                yv->SetName(var);

                // Assemble the curve.
                int nPts = mesh->second[0];
                vtkRectilinearGrid *rg = vtkVisItUtility::Create1DRGrid(nPts, VTK_FLOAT);
                vtkFloatArray *xc = vtkFloatArray::SafeDownCast(rg->GetXCoordinates());
                bool validCoords = coordvals && (sz == static_cast<size_t>(nPts));
                for (int j = 0 ; j < nPts ; ++j)
                    xc->SetValue(j, validCoords ? coordvals[j] : (float)j);

                rg->GetPointData()->SetScalars(yv);
                retval = rg;

                delete [] coordvals;
                yv->Delete();
            }
            else
            {
                // Static curves

                int nPts = mesh->second[0];
                vtkFloatArray *yv = vtkFloatArray::New();
                yv->SetNumberOfTuples(nPts); 
                if(fileObject->ReadVariableIntoAsFloat(var, (float*)yv->GetVoidPointer(0)))
                {
                    vtkRectilinearGrid *rg = vtkVisItUtility::Create1DRGrid(
                                             nPts, VTK_FLOAT);
                    vtkFloatArray *xc = vtkFloatArray::SafeDownCast(
                                            rg->GetXCoordinates());
                    yv->SetName(var);
                    for (int j = 0 ; j < nPts ; j++)
                        xc->SetValue(j, (float)j); 

                    rg->GetPointData()->SetScalars(yv);
                    retval = rg;
                }

                yv->Delete();
            }
        }
        else
        {
            vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New(); 

            // Get the dim starts and counts.
            intVector dimStarts, dimCounts;
            ReturnDimStartsAndCounts(timeState, mesh->second, dimStarts, dimCounts);

            // Get the indices of the spatial dimensions
            int spatialDimIndices[3] = {0,0,0}, nSpatialDims = 0;
            ReturnSpatialDimensionIndices(mesh->second, spatialDimIndices, nSpatialDims);
             
            //
            // Populate the coordinates.  Put in 3D points with z=0 if
            // the mesh is 2D.
            //
            const char *xyzname[] = {"X", "Y", "Z"};
            int            dims[3] = {0,0,0};
            vtkFloatArray *coords[3] = {0,0,0};
            for (int i = 0 ; i < 3 ; i++)
            {
                // Default number of components for an array is 1.
                coords[i] = vtkFloatArray::New();

                if (i < nSpatialDims)
                {
                    // Get the name of this dimension and try and read an array by the same 
                    // name, if one exists and is the same size as the dimension. We'll use 
                    // it for the mesh's coordinate values.
                    float *coordvals = 0;
                    char dimName[NC_MAX_NAME+1];
                    StringIntVectorMap::const_iterator meshDims = meshNameToNCDimensions.find(bare_var);
                    if(meshDims != meshNameToNCDimensions.end())
                    {
                        int status;
                        size_t sz;
                        debug4 << mName << "Looking for " << xyzname[i]
                               << " coordinate dimension name" << endl;
                        if((status = nc_inq_dim(fileObject->GetFileHandle(), 
                            meshDims->second[spatialDimIndices[i]], dimName, &sz)) == NC_NOERR)
                        {
                            TypeEnum dvart = NO_TYPE;
                            int dvarndims = 0;
                            int *dvardims = 0;
                            debug4 << mName << "Looking for " << dimName << " array as "
                                   << xyzname[i] << " coordinate" << endl;
                            if(fileObject->InqVariable(dimName, &dvart, &dvarndims, &dvardims))
                            {
                                if(dvarndims == 1 && (size_t)dvardims[0] == sz)
                                    coordvals = ReadArray(dimName);
                                delete [] dvardims;
                            }
                        }
                    }

                    // Populate the coordinates
                     dims[i] = dimCounts[spatialDimIndices[i]];
                    if (centering == AVT_NODECENT)
                        coords[i]->SetNumberOfTuples(dims[i]);
                    else
                        coords[i]->SetNumberOfTuples(dims[i]+1);
                    if(coordvals != 0)
                    {
                        debug4 << mName << "Using " << dimName << " array as "
                                   << xyzname[i] << " coordinate" << endl;
                        int start = dimStarts[spatialDimIndices[i]];
                        if (centering == AVT_NODECENT)
                        {
                            for (int j = 0 ; j < dims[i] ; j++)
                                coords[i]->SetComponent(j, 0, coordvals[start + j]);
                        }
                        else // mesh coords for zone-centered variants
                        {
                            float width_of_first_zone =
                                coordvals[start+1] - coordvals[start];
                            float half_step_before_first =
                                coordvals[start] - width_of_first_zone / 2.0;
                            coords[i]->SetComponent(0, 0, half_step_before_first);

                            // Coords for zone-centered variant are midway between
                            // coords for the normal case
                            for (int j = 0 ; j < dims[i]-1 ; j++)
                            {
                                float midway =
                                    (coordvals[start+j] + coordvals[start+j+1]) / 2.0;
                                coords[i]->SetComponent(j+1, 0, midway);
                            }

                            float width_of_last_zone =
                                coordvals[start+dims[i]-1] - coordvals[start+dims[i]-2];
                            float half_step_after_last =
                                coordvals[start+dims[i]-1] + width_of_last_zone / 2.0;
                            coords[i]->SetComponent(dims[i], 0, half_step_after_last);

                            // Don't forget to adjust dims sizes up by 1
                            dims[i] = dimCounts[spatialDimIndices[i]]+1;
                        }
                    }
                    else
                    {
                        debug4 << mName << "Using indices for " << xyzname[i] << " coordinate" << endl;
                        for (int j = 0 ; j < dims[i] ; j++)
                            coords[i]->SetComponent(j, 0, dimStarts[spatialDimIndices[i]] + j);
                    }

                    if(coordvals != 0)
                        delete [] coordvals;
                }
                else
                {
                    dims[i] = 1;
                    coords[i]->SetNumberOfTuples(1);
                    coords[i]->SetComponent(0, 0, 0.);
                }
            }
            rgrid->SetDimensions(dims);
            rgrid->SetXCoordinates(coords[0]);
            coords[0]->Delete();
            rgrid->SetYCoordinates(coords[1]);
            coords[1]->Delete();
            rgrid->SetZCoordinates(coords[2]);
            coords[2]->Delete();

            vtkIntArray *arr = vtkIntArray::New();
            arr->SetNumberOfTuples(3);
            arr->SetValue(0, dimStarts[0]);
            arr->SetValue(1, dimStarts[1]);
            arr->SetValue(2, (dimStarts.size() > 2 ? dimStarts[2] : 0));
            arr->SetName("base_index");
            rgrid->GetFieldData()->AddArray(arr);
            arr->Delete();

            retval = rgrid;
        }
    }
    else
    {
        EXCEPTION1(InvalidVariableException, bare_var);
    }

    return retval;
}

// ****************************************************************************
// Method: avtBasicNETCDFReader::GetVar
//
// Purpose: 
//   Returns the data for the specified variable.
//
// Arguments:
//   timeState : The time state to read.
//   var       : The name of the variable to read.
//
// Returns:    The data or 0.
//
// Note:       
//
// Programmer: Brad Whitlock
// Creation:   Thu Aug 18 18:06:49 PST 2005
//
// Modifications:
//
//   Mark C. Miller, Tue Aug 15 15:28:11 PDT 2006
//   Added logic to support on-the-fly domain decomposition. Added macro to
//   do a partial read
//
//   Mark C. Miller, Wed Aug 16 14:45:22 PDT 2006
//   Fixed reversal of coords between VisIt and netcdf
//   
//   Jeremy Meredith, Thu Mar 22 15:05:02 EDT 2007
//   Since we're calling InqVariable here, we need to take the same steps
//   as we do when reading the mesh dims to remove trailing 1's.  Also,
//   print to debug4 the valid dimensions array instead of the raw one.
//
//   Jeremy Meredith, Wed Apr  4 14:01:40 EDT 2007
//   Adjusted to account for an occasional "4-dimensional" mesh.
//
//   Brad Whitlock, Thu Oct 29 15:23:32 PDT 2009
//   I rewrote the code so it can read variables that have a time dimension.
//
//   Kathleen Bonnell, Wed Sep 8 20:34:11 MST 2010 
//   Compare minfo to varToDimensionSizes.end, not meshNameToDimensionsSizes.
//
//   Brad Whitlock, Fri Jan  6 13:35:00 PST 2012
//   Apply scaling and offset attributes to transform data.
//
//   Mark C. Miller, Wed Nov 20 14:47:30 PST 2013
//   Handle zone centered variants of variables too.
// ****************************************************************************

#define READVAR(VTKTYPE) \
        {\
            VTKTYPE *arr = VTKTYPE::New();\
            arr->SetNumberOfComponents(1);\
            arr->SetNumberOfTuples(nValues);\
            debug4 << "Allocated a " << \
                    #VTKTYPE \
                   << " of " << nValues << " elements" << endl; \
            int *rdimStarts = new int[dimStarts.size()]; \
            int *rdimCounts = new int[dimCounts.size()]; \
            for (size_t kk = 0; kk < dimStarts.size(); kk++)\
            {\
                rdimStarts[ndims-kk-1] = dimStarts[kk];\
                rdimCounts[ndims-kk-1] = dimCounts[kk];\
            }\
            if(fileObject->ReadVariableInto(bare_var, t, rdimStarts, rdimCounts,\
                                            arr->GetVoidPointer(0)))\
                retval = arr;\
            else\
                arr->Delete();\
            delete [] rdimStarts;\
            delete [] rdimCounts;\
        }


vtkDataArray *
avtBasicNETCDFReader::GetVar(int timeState, const char *var)
{
    const char *mName = "avtBasicNETCDFReader::GetVar: ";
    debug4 << mName << "var=" << var << endl;
    char const *bare_var = var;
    if (!strncmp(var, "as_zonal/", sizeof("as_zonal/")-1))
        bare_var = &var[sizeof("as_zonal/")-1];

    vtkDataArray *retval = 0;

    // Determine the storage type for the variable.
    TypeEnum t = NO_TYPE;
    int ndims = 0;
    int *dims = 0;
    if(!fileObject->InqVariable(bare_var, &t, &ndims, &dims))
    {
        EXCEPTION1(InvalidVariableException, var);
    }
    delete [] dims;

    // Look up the mesh dimensions for the mesh.
    StringIntVectorMap::const_iterator minfo = varToDimensionsSizes.find(bare_var);
    if(minfo != varToDimensionsSizes.end())
    {
        // Get the dim starts and counts.
        intVector dimStarts, dimCounts;
        ReturnDimStartsAndCounts(timeState, minfo->second, dimStarts, dimCounts);

        // Compute how many values that makes.
        unsigned long nValues = 1;
        for(size_t i = 0; i < dimCounts.size(); ++i)
            nValues *= (unsigned long)dimCounts[i];

        debug4 << mName << "dimStarts = {";
        for(size_t i = 0; i < dimStarts.size(); ++i) {
            debug4 << dimStarts[i] << ", ";
        }
        debug4 << "}\n";
        debug4 << mName << "dimCounts = {";
        for(size_t i = 0; i < dimCounts.size(); ++i) {
            debug4 << dimCounts[i] << ", ";
        }
        debug4 << "}\n";

        if(t == CHARARRAY_TYPE || t == UCHARARRAY_TYPE)
        {
            READVAR(vtkUnsignedCharArray)

            vtkDataArray *scaled = ApplyScaleAndOffset(var, retval);
            retval->Delete();
            retval = scaled;
        }
        else if(t == SHORTARRAY_TYPE)
        {
            READVAR(vtkShortArray)

            vtkDataArray *scaled = ApplyScaleAndOffset(var, retval);
            retval->Delete();
            retval = scaled;
        }
        else if(t == INTEGERARRAY_TYPE)
            READVAR(vtkIntArray)
        else if(t == LONGARRAY_TYPE)
            READVAR(vtkLongArray)
        else if(t == FLOATARRAY_TYPE)
            READVAR(vtkFloatArray)
        else if(t == DOUBLEARRAY_TYPE)
            READVAR(vtkDoubleArray)
    }
    else
    {
        EXCEPTION1(InvalidVariableException, var);
    }

    return retval;
}

