/*****************************************************************************
*
* Copyright (c) 2000 - 2018, 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.
*
*****************************************************************************/

// ************************************************************************* //
//                         avtPICS_TesterFileFormat.C                        //
// ************************************************************************* //

#include <avtPICS_TesterFileFormat.h>

#include <string>

#include <vtkDoubleArray.h>
#include <vtkPoints.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkFieldData.h>

#include <avtDatabase.h>
#include <avtDatabaseMetaData.h>
#include <avtIntervalTree.h>
#include <avtStructuredDomainBoundaries.h>
#include <avtVariableCache.h>

#include <avtParallel.h>

#include <DBOptionsAttributes.h>
#include <Expression.h>

#include <DebugStream.h>
#include <InvalidVariableException.h>

// ****************************************************************************
//  Method: PrintVec
//
//  Purpose:
//      Convenience function.  Print a vector. 
//
//  Programmer: Dave Bremer
//  Creation:   Fri Jun 13 15:54:11 PDT 2008
//
// ****************************************************************************
template <class T>
std::string Vec2String(std::string name, T *vec, int numelems) {
  std::string s = name + " = [" ;
  int elem = 0;
  char buf[32];
  
  while (elem < numelems ) {
    float value = vec[elem];
    SNPRINTF(buf,31,"%f",value); 
    s += buf ;
    if (elem == numelems - 1) {
      s+= "]"; 
    } else {
      s+= ", "; 
    }
    elem ++; 
  }
  return s;
} 


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat constructor
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

avtPICS_TesterFileFormat::avtPICS_TesterFileFormat(const char *filename)
    : avtMTMDFileFormat(filename)
{
    flowType = STANDARD;
    rank = 3;
    isRectilinear = true;

    ReadHeader(filename);
}

void
avtPICS_TesterFileFormat::ReadHeader(const char *filename)
{
    times.clear();
    cycles.clear();
    
    numBlocks[0].clear();
    numBlocks[1].clear();
    numBlocks[2].clear();
    
    numCells[0].clear();
    numCells[1].clear();
    numCells[2].clear();

    char line[1024];
    ifstream ifile(filename);

    while (! ifile.eof())
    {
      ifile.getline(line, 1024);

      if (line[0] == '#' || strlen(line) == 0)
        continue;
      else if (strcmp(line, "3D") == 0)
      {
        rank = 3;
      }
      else if (strcmp(line, "2D") == 0)
      {
       rank = 2;
      }
      else if (strcmp(line, "rectilinear") == 0)
      {
        flowType = STANDARD;
        isRectilinear = true;
      }
      else if (strcmp(line, "unstructured") == 0)
      {
        flowType = STANDARD;
        isRectilinear = false;
      }
      else if (strcmp(line, "DOUBLE_GYRE") == 0)
      {
        flowType = DOUBLE_GYRE;
        rank = 3;
        isRectilinear = true;
      }
      else if (strcmp(line, "ABC_FLOW_STEADY_STATE") == 0)
      {
        flowType = ABC_FLOW_STEADY_STATE;
        rank = 3;
        isRectilinear = true;
      }
      else if (strcmp(line, "ABC_FLOW_APERIODIC") == 0)
      {
        flowType = ABC_FLOW_APERIODIC;
        rank = 3;
        isRectilinear = true;
      }
      
      else if( flowType == STANDARD )
      {
        double time;
        int nb;
        int nc;
        double velx, vely, velz;
        double magnitude;
        bool validScan = false;
        
        global_extents[0] = 1.0;
        global_extents[1] = 1.0;

        if (rank == 2)
          global_extents[2] = 0.0;
        else
          global_extents[2] = 1.0;
          
        if (rank == 3)
        {
         validScan = (sscanf(line, "%lf %d %d %lf %lf %lf %lf",
                             &time, &nb, &nc,
                             &velx, &vely, &velz, &magnitude) == 7);
        }
        else
        {
         validScan = (sscanf(line, "%lf %d %d %lf %lf %lf",
                             &time, &nb, &nc,
                             &velx, &vely, &magnitude) == 6);
         velz = 0.;
        }
 
        if (validScan)
        {
         times.push_back(time);
         cycles.push_back(cycles.size());
         
         numBlocks[0].push_back(nb);
         numBlocks[1].push_back(nb);
         
         if (rank == 2)
           numBlocks[2].push_back(1);
         else
           numBlocks[2].push_back(nb);
         
         numCells[0].push_back(nc);
         numCells[1].push_back(nc);
         
         if (rank == 2)
           numCells[2].push_back(1);
         else
           numCells[2].push_back(nc);
         
         vels.push_back(velx);
         vels.push_back(vely);
         vels.push_back(velz);
         
         magnitudes.push_back(magnitude);
        }
        else
          flowType = UNKNOWN;
      }
      else if( flowType == DOUBLE_GYRE )
      {
        int nTimes, nb, nx, ny, nz;
        double startTime, endTime, simulationTime, dt, bx, by, bz;
        
        bool validScan =
          (sscanf(line, "%d %lf %lf %d %d %d %d %lf %lf %lf %lf %lf %lf",
                  &nTimes, &startTime, &simulationTime,
                  &nb, &nx, &ny, &nz,
                  &bx, &by, &bz,
                  &dg_A, &dg_epsilon, &dg_period) == 13);
 
        if( validScan )
        {
          times.clear();
          cycles.clear();

          numBlocks[0].clear();
          numBlocks[1].clear();
          numBlocks[2].clear();
          
          numCells[0].clear();
          numCells[1].clear();
          numCells[2].clear();

          if (nz == 1)
            rank = 2;
          else
            rank = 3;
              
          global_extents[0] = bx;
          global_extents[1] = by;
          global_extents[2] = bz;

          if( simulationTime > 0 )
            endTime = startTime + simulationTime;
          else
          {
            endTime = startTime;
            startTime = endTime + simulationTime;
          }

          dt = (endTime - startTime) / (double) nTimes;

          for( int i=0; i<nTimes+1; ++i )
          {
            times.push_back( startTime + (double) i * dt );
            cycles.push_back(i);
            
            numBlocks[0].push_back(nb);
            numBlocks[1].push_back(nb);

            if (rank == 2)
              numBlocks[2].push_back(1);
            else
              numBlocks[2].push_back(nb);
         
            numCells[0].push_back(nx);
            numCells[1].push_back(ny);

            if (rank == 2)
              numCells[2].push_back(1);
            else
              numCells[2].push_back(nz);
          }
        }
        else
          flowType = UNKNOWN;
      }
      else if( flowType == ABC_FLOW_STEADY_STATE )
      {
        int nb, nx, ny, nz;

        bool validScan = (sscanf(line, "%d %d %d %d",
                                 &nb, &nx, &ny, &nz ) == 4);
 
        if( validScan )
        {
          times.clear();
          cycles.clear();

          numBlocks[0].clear();
          numBlocks[1].clear();
          numBlocks[2].clear();
          
          numCells[0].clear();
          numCells[1].clear();
          numCells[2].clear();

          global_extents[0] = 2.0 * M_PI;
          global_extents[1] = 2.0 * M_PI;
          global_extents[2] = 2.0 * M_PI;
          
          times.push_back(0);       
          cycles.push_back(0);
            
          numBlocks[0].push_back(nb);
          numBlocks[1].push_back(nb);
          numBlocks[2].push_back(nb);
         
          numCells[0].push_back(nx);
          numCells[1].push_back(ny);
          numCells[2].push_back(nz);
        }
        else
          flowType = UNKNOWN;
      }
      else if( flowType == ABC_FLOW_APERIODIC )
      {
        int nTimes, nb, nx, ny, nz;
        double startTime, endTime, simulationTime, dt;
        
        bool validScan =
          (sscanf(line, "%d %lf %lf %d %d %d %d %lf %lf %lf %d %d %d",
                  &nTimes, &startTime, &simulationTime,
                  &nb, &nx, &ny, &nz,
                  &abc_c0, &abc_c1, &abc_c2,
                  &abc_signalA, &abc_signalB, &abc_signalC) == 13);
 
        if( validScan )
        {
          times.clear();
          cycles.clear();

          numBlocks[0].clear();
          numBlocks[1].clear();
          numBlocks[2].clear();
          
          numCells[0].clear();
          numCells[1].clear();
          numCells[2].clear();

          global_extents[0] = 2.0 * M_PI;
          global_extents[1] = 2.0 * M_PI;
          global_extents[2] = 2.0 * M_PI;
          
          if( simulationTime > 0 )
            endTime = startTime + simulationTime;
          else
          {
            endTime = startTime;
            startTime = endTime + simulationTime;
          }

          dt = (endTime - startTime) / (double) nTimes;

          for( int i=0; i<nTimes+1; ++i )
          {
            times.push_back( startTime + (double) i * dt );
            cycles.push_back(i);
            
            numBlocks[0].push_back(1);
            numBlocks[1].push_back(1);
            numBlocks[2].push_back(1);
            
            numCells[0].push_back(nx);
            numCells[1].push_back(ny);
            numCells[2].push_back(nz);
          }
        }
        else
          flowType = UNKNOWN;
      }
    }
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

int
avtPICS_TesterFileFormat::GetNTimesteps(void)
{
    return times.size();
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::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: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

void
avtPICS_TesterFileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::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: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

void
avtPICS_TesterFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timestate)
{
    std::string meshname = "mesh";
    avtMeshType meshtype = AVT_RECTILINEAR_MESH;
    if (! isRectilinear)
        meshtype = AVT_UNSTRUCTURED_MESH;

    double extents[6];
    extents[0] = 0;
    extents[1] = global_extents[0];
    extents[2] = 0;
    extents[3] = global_extents[1];
    extents[4] = 0;
    extents[5] = global_extents[2];

    int nblocks = 1;

    for( int i=0; i<rank; ++i )
      nblocks *= numBlocks[i][timestate];

    int block_origin = 0;
    int spatial_dimension = rank;
    int topological_dimension = rank;

    int bounds[3] = {numCells[0][timestate]+1,
                     numCells[1][timestate]+1,
                     numCells[2][timestate]+1};

    if (rank == 2)
      bounds[2] = 1;

    avtMeshMetaData *mesh = new avtMeshMetaData;
    mesh->name = meshname;
    mesh->meshType = meshtype;
    mesh->numBlocks = nblocks;
    mesh->blockOrigin = block_origin;
    mesh->cellOrigin = 0;
    mesh->spatialDimension = spatial_dimension;
    mesh->topologicalDimension = topological_dimension;
    mesh->blockTitle = "blocks";
    mesh->blockPieceName = "block";
    mesh->SetBounds(bounds);
    mesh->hasLogicalBounds = true;
    mesh->SetExtents(extents);
    mesh->hasSpatialExtents = true;
    mesh->containsGhostZones = AVT_NO_GHOSTS;

    md->Add(mesh);


    std::string varname = "velocity";
    int vector_dim = spatial_dimension;
    avtCentering cent = AVT_NODECENT;
    AddVectorVarToMetaData(md, varname, meshname, cent, vector_dim);

    md->SetTimes(times);
    md->SetTimesAreAccurate(true);
    md->SetCycles(cycles);
    md->SetCyclesAreAccurate(true);

    md->SetTemporalExtents(times[0], times[times.size()-1]);
    md->SetHasTemporalExtents(true);

    // Find logical domain boundaries
    if (!avtDatabase::OnlyServeUpMetaData() && nblocks > 1)
    {
        avtRectilinearDomainBoundaries *rdb =
          new avtRectilinearDomainBoundaries(true);

        rdb->SetNumDomains(nblocks);
        int bbox[6];

        for (int domain = 0; domain < (size_t)nblocks ; ++domain)
        {
            int xOff = domain % numBlocks[0][timestate];
            int yOff = (domain/numBlocks[0][timestate]) % numBlocks[1][timestate];
            int zOff = domain/(numBlocks[0][timestate]*numBlocks[1][timestate]);

            bbox[0] = (xOff    ) * numCells[0][timestate];
            bbox[1] = (xOff + 1) * numCells[0][timestate];

            bbox[2] = (yOff    ) * numCells[1][timestate];
            bbox[3] = (yOff + 1) * numCells[1][timestate];

            // VisIt expects the 2d case to have flat logical z extent (0,0).
            if(rank == 2)
            {
                bbox[4] = 0;
                bbox[5] = 0;
            }
            else // if(rank == 3)
            {
              bbox[4] = (zOff    ) * numCells[2][timestate];
              bbox[5] = (zOff + 1) * numCells[2][timestate];
            }

            rdb->SetIndicesForRectGrid(domain, bbox);
        }

        rdb->CalculateBoundaries();

        void_ref_ptr vr =
          void_ref_ptr(rdb, avtRectilinearDomainBoundaries::Destruct);

        cache->CacheVoidRef("any_mesh",
                            AUXILIARY_DATA_DOMAIN_BOUNDARY_INFORMATION,
                            timestate, -1, vr);
    }
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::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.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

void
RotatePoint(const double *pt3, double *pt2)
{
    double mid[3] = { 0.5, 0.5, 0.5 };

    double pt[3];
    pt[0] = pt3[0] - mid[0];
    pt[1] = pt3[1] - mid[1];
    pt[2] = pt3[2] - mid[2];

    pt2[0] = pt[0]*0.707-pt[1]*0.707; // x*cos\theta-y*sin\theta, theta=45
    pt2[1] = pt[0]*0.707+pt[1]*0.707; // x*sin\theta+y*cos\theta, theta=45
    pt2[2] = pt[2];

    pt2[0] += mid[0];
    pt2[1] += mid[1];
    pt2[2] += mid[2];
}

vtkDataSet *
avtPICS_TesterFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
    int nblocks = 1;

    for( int i=0; i<rank; ++i )
      nblocks *= numBlocks[i][timestate];

    if (domain >= nblocks)
    {
        EXCEPTION1(VisItException, "Invalid mesh requested!");
    }

    double xSizePerBlock = global_extents[0]/numBlocks[0][timestate];
    double ySizePerBlock = global_extents[1]/numBlocks[1][timestate];
    double zSizePerBlock = global_extents[2]/numBlocks[2][timestate];

    int xOff = domain % numBlocks[0][timestate];
    int yOff = (domain/numBlocks[0][timestate]) % numBlocks[1][timestate];
    int zOff = domain/(numBlocks[0][timestate]*numBlocks[1][timestate]);

    int dims[3] = { numCells[0][timestate]+1,
                    numCells[1][timestate]+1,
                    numCells[2][timestate]+1 };
    if (rank == 2)
      dims[2] = 1;

    if (isRectilinear)
    {
        vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
    
        rgrid->SetDimensions(dims);
    
        vtkDoubleArray *x = vtkDoubleArray::New();
        x->SetNumberOfTuples(numCells[0][timestate]+1);
        for (int i = 0 ; i < numCells[0][timestate]+1 ; i++)
            x->SetTuple1(i, xSizePerBlock*xOff + i*xSizePerBlock/numCells[0][timestate]);
        rgrid->SetXCoordinates(x);
        x->Delete();
    
        vtkDoubleArray *y = vtkDoubleArray::New();
        y->SetNumberOfTuples(numCells[1][timestate]+1);
        for (int i = 0 ; i < numCells[1][timestate]+1 ; i++)
            y->SetTuple1(i, ySizePerBlock*yOff + i*ySizePerBlock/numCells[1][timestate]);
        rgrid->SetYCoordinates(y);
        y->Delete();
    
        if (rank == 3)
        {
            vtkDoubleArray *z = vtkDoubleArray::New();
            z->SetNumberOfTuples(numCells[2][timestate]+1);
            for (int i = 0 ; i < numCells[2][timestate]+1 ; i++)
                z->SetTuple1(i, zSizePerBlock*zOff + i*zSizePerBlock/numCells[2][timestate]);
            rgrid->SetZCoordinates(z);
            z->Delete();
        } //if (rank == 2)
        else
        {
            vtkDoubleArray *z = vtkDoubleArray::New();
            z->SetNumberOfTuples(1);
            z->SetTuple1(0, 0.0);
            rgrid->SetZCoordinates(z);
            z->Delete();
        }

        // Sneak in periodic boundaries for the avtIVPsolver.
        if( flowType == ABC_FLOW_STEADY_STATE ||
            flowType == ABC_FLOW_APERIODIC )
        {
            vtkDoubleArray *bounds = vtkDoubleArray::New();
            bounds->SetName("Periodic Boundaries");

            // Set the number of components before setting the number of tuples
            // for proper memory allocation.
            bounds->SetNumberOfComponents( 1 );
            bounds->SetNumberOfTuples( 3 );

            bounds->SetTuple1(0, global_extents[0] );
            bounds->SetTuple1(1, global_extents[1] );
            bounds->SetTuple1(2, global_extents[2] );
  
            rgrid->GetFieldData()->AddArray(bounds);

            bounds->Delete();
        }
    
        return rgrid;
    }
    else
    {
        vtkPoints *pts = vtkPoints::New();
        //pts->SetNumberOfComponents(3);
        int npts = dims[0]*dims[1]*dims[2];
        pts->SetNumberOfPoints(npts);

        int idx = 0;
        for (int i = 0 ; i < dims[0] ; i++)
            for (int j = 0 ; j < dims[1] ; j++)
                for (int k = 0 ; k < dims[2] ; k++)
                {
                    double pt[3];
                    pt[0] = xSizePerBlock*xOff +
                      i*xSizePerBlock/numCells[0][timestate];
                    pt[1] = ySizePerBlock*yOff +
                      j*ySizePerBlock/numCells[1][timestate];
                    pt[2] = zSizePerBlock*zOff +
                      k*zSizePerBlock/numCells[2][timestate];
  
                    double pt2[3];
                    RotatePoint(pt, pt2);
                    pts->SetPoint(idx++, pt2);
                }


        vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
        sgrid->SetDimensions(dims);
        sgrid->SetPoints(pts);
        pts->Delete();

        return sgrid;
    }
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, 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.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

vtkDataArray *
avtPICS_TesterFileFormat::GetVar(int timestate, int domain, const char *varname)
{
    return NULL;
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkDoubleArray, 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.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

vtkDataArray *
avtPICS_TesterFileFormat::GetVectorVar(int timestate, int domain, const char *varname)
{
    double t = times[timestate];

    double pt[3], vec[3];

    int dims[3] = { numCells[0][timestate]+1,
                    numCells[1][timestate]+1,
                    numCells[2][timestate]+1 };
    if (rank == 2)
        dims[2] = 1;

    int ntuples = dims[0]*dims[1]*dims[2];
    vtkDoubleArray *rv = vtkDoubleArray::New();
    rv->SetNumberOfComponents(3);
    rv->SetNumberOfTuples(ntuples);

    if( flowType == STANDARD )
    {
        srand((timestate+1)*(domain+1));
        for (int i = 0 ; i < ntuples ; i++)
        {
            for (int j = 0 ; j < 3 ; j++)
            {
                if( j < rank )
                {
                    double r = (rand()%1000)/1000.0;
                    vec[j] = vels[3*timestate+j]+magnitudes[timestate]*r;
                }
                else
                  vec[j] = 0.0;
            }
            rv->SetTuple(i, vec);
        }
    }
    else if( flowType == DOUBLE_GYRE )
    {
        vtkDataSet * rectGrid = GetMesh(timestate, domain, "mesh");

        double A = dg_A;
        double epsilon = dg_epsilon;
        double omega = 2.0 * M_PI / dg_period;
            
        for (int i = 0 ; i < ntuples ; i++)
        {
            rectGrid->GetPoint( i, pt );

            // Test code for a double gyre.
            double xi = pt[0];
            double yi = pt[1];
            
            double at = epsilon * sin( omega * t );
            double bt = 1.0  - 2.0 * at;
            
            double fxt = (at * xi + bt) * xi;
            double dfx = (2.0 * at * xi + bt);
            
            vec[0] = -M_PI * A * sin(M_PI * fxt) * cos(M_PI * yi);
            vec[1] =  M_PI * A * cos(M_PI * fxt) * sin(M_PI * yi) * dfx;
            vec[2] = 0;

            rv->SetTuple(i, vec);
        }

        rectGrid->Delete();
    }
    else if( flowType == ABC_FLOW_STEADY_STATE )
    {
        vtkDataSet * rectGrid = GetMesh(timestate, domain, "mesh");

        double A = sqrt(3.0);
        double B = sqrt(2.0);
        double C = 1.0;
 
        for (int i = 0 ; i < ntuples ; i++)
        {
            rectGrid->GetPoint( i, pt );

            // Test code for the ABC.
            double xi = pt[0];
            double yi = pt[1];
            double zi = pt[2];
            
            vec[0] = A * sin(zi) + C * cos(yi);
            vec[1] = B * sin(xi) + A * cos(zi);
            vec[2] = C * sin(yi) + B * cos(xi);

            rv->SetTuple(i, vec);
        }

        rectGrid->Delete();
    }
    else if( flowType == ABC_FLOW_APERIODIC )
    {
        vtkDataSet * rectGrid = GetMesh(timestate, domain, "mesh");

        double A = sqrt(3.0);
        double B = sqrt(2.0);
        double C = 1.0;
        
        double c0 = abc_c0;
        double c1 = abc_c1;
        double c2 = abc_c2;

        for (int i = 0 ; i < ntuples ; i++)
        {
            rectGrid->GetPoint( i, pt );

            // Test code for the ABC.
            double xi = pt[0];
            double yi = pt[1];
            double zi = pt[2];

            double signalA =
              (abc_signalA ? A*c0*tanh(c1*t)*sin((c2*t)*(c2*t)) : 0);
            double signalB = 
              (abc_signalB ? B*c0*tanh(c1*t)*cos((c2*t)*(c2*t)) : 0);
            double signalC = 
              (abc_signalC ? C*c0*tanh(c1*t)*sin((c2*t)*(c2*t)) : 0);

            vec[0] = (A+signalA) * sin(zi) + (C+signalC) * cos(yi);
            vec[1] = (B+signalB) * sin(xi) + (A+signalA) * cos(zi);
            vec[2] = (C+signalC) * sin(yi) + (B+signalB) * cos(xi);

            rv->SetTuple(i, vec);
        }

        rectGrid->Delete();
    }

    return rv;
}


// ****************************************************************************
//  Method: avtPICS_TesterFileFormat::GetAuxiliaryData
//
//  Purpose:
//      Gets the auxiliary data from a Silo file.
//
//  Arguments:
//      var        The variable of interest.
//      domain     The domain of interest.
//      type       The type of auxiliary data.
//      <unnamed>  The arguments for that -- not used for any PICS types.
//      df         The interval tree destructor function.
//
//  Returns:    The auxiliary data.  Throws an exception if this is not a
//              supported data type.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Tue Mar 6 07:45:19 PDT 2012
//
// ****************************************************************************

void *
avtPICS_TesterFileFormat::GetAuxiliaryData(const char *var, int ts, int dom,
                                           const char * type, void *,
                                           DestructorFunction &df)
{
    if (strcmp(type, AUXILIARY_DATA_SPATIAL_EXTENTS) == 0)
    {
        int nblocks = 1;

        for( int i=0; i<rank; ++i )
          nblocks *= numBlocks[i][ts];

        int dimension = rank;
        avtIntervalTree *itree = new avtIntervalTree(nblocks, dimension);

        double xSizePerBlock = global_extents[0]/numBlocks[0][ts];
        double ySizePerBlock = global_extents[1]/numBlocks[1][ts];
        double zSizePerBlock = global_extents[2]/numBlocks[2][ts];

        for (int domain = 0 ; domain < nblocks ; domain++)
        {
            int xOff = domain % numBlocks[0][ts];
            int yOff = (domain/numBlocks[0][ts]) % numBlocks[1][ts];
            int zOff = domain/(numBlocks[0][ts]*numBlocks[1][ts]);

            double extents[6];
            extents[0] = xOff*xSizePerBlock;
            extents[1] = (xOff+1)*xSizePerBlock;
            extents[2] = yOff*ySizePerBlock;
            extents[3] = (yOff+1)*ySizePerBlock;

            if (rank == 3)
            {
                extents[4] = zOff*zSizePerBlock;
                extents[5] = (zOff+1)*zSizePerBlock;
            }
            else
            {
                extents[4] = 0.0;
                extents[5] = 0.0;
            }
            if (! isRectilinear)
            {
                int i;
                double b[6] = { 1.0e6, -.10e6, 1.0e6, -1.0e6, 1.0e6, -1.0e6 };

                if (rank == 2)
                   b[4] = b[5] = 0.0;

                for (i = 0 ; i < 8 ; i++)
                {
                    double pt[3];
                    pt[0] = (i%2 ? extents[0] : extents[1]);
                    pt[1] = (((i/2)%2) ? extents[2] : extents[3]);
                    pt[2] = (i/2 ? extents[4] : extents[5]);
                    double pt2[3];
                    RotatePoint(pt, pt2);
                    b[0] = (b[0] < pt2[0] ? b[0] : pt2[0]);
                    b[1] = (b[1] > pt2[0] ? b[1] : pt2[0]);
                    b[2] = (b[2] < pt2[1] ? b[2] : pt2[1]);
                    b[3] = (b[3] > pt2[1] ? b[3] : pt2[1]);
                    b[4] = (b[4] < pt2[2] ? b[4] : pt2[2]);
                    b[5] = (b[5] > pt2[2] ? b[5] : pt2[2]);
                }
                for (i = 0 ; i < 6 ; i++)
                    extents[i] = b[i];
            }
            itree->AddElement(domain, extents);
        }
        itree->Calculate(true);

        df = avtIntervalTree::Destruct;

        return ((void *) itree);
    }

    return NULL;
}
