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

// ************************************************************************* //
//                            avtCarpetN5FileFormat.C                           //
// ************************************************************************* //

#include <avtCarpetN5FileFormat.h>

#include <string>

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

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

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

#include <InvalidVariableException.h>



// ****************************************************************************
//  Method: avtCarpetN5FileFormat constructor
//
//  Programmer: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

avtCarpetN5FileFormat::avtCarpetN5FileFormat(const char *filename)
    : File(filename), avtMTMDFileFormat(filename), nlevels(0)
{
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

int
avtCarpetN5FileFormat::GetNTimesteps(void)
{
    return File.GetCycles().size();
}


void avtCarpetN5FileFormat::GetCycles(std::vector<int>& c)
{
   c = File.GetCycles();
}



void avtCarpetN5FileFormat::GetTimes(std::vector<double>& t)
{
   t = File.GetTimes();
}


// ****************************************************************************
//  Method: avtCarpetN5FileFormat::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: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

void
avtCarpetN5FileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avtCarpetN5FileFormat::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: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

void
avtCarpetN5FileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
    md->SetTimes(File.GetTimes());
    md->SetTimesAreAccurate(true);
    md->SetCycles(File.GetCycles());
    md->SetCyclesAreAccurate(true);

    
    // loop over all meshes at current timestate and timelevel=0 contained in file
    Nirvana::metadata   select_by_identifier;
    Nirvana::metadata   select_by_value;
    select_by_identifier << attribute<int>("cycle", 0);
    select_by_identifier << attribute<int>("timelevel", 0);
    select_by_identifier << attribute<string>("meshname", "");
    select_by_value      << attribute<int>   ("cycle", File.GetCycles()[timeState]);
    select_by_value      << attribute<int>   ("timelevel", 0);
    for (CarpetN5::Iterator iMeshes = File.Begin(select_by_identifier, select_by_value); !iMeshes.Done(); ++iMeshes)
    {
       Nirvana::metadata   select_by_identifier;
       Nirvana::metadata   select_by_value;
       
       select_by_identifier << attribute<int>("cycle", 0);
       select_by_identifier << attribute<int>("timelevel", 0);
       select_by_identifier << attribute<string>("meshname", "");
       select_by_identifier << attribute<int>("map", 0);
       select_by_identifier << attribute<int>("reflevel", 0);
       select_by_identifier << attribute<int>("comp", 0);
       
       select_by_value      << (*iMeshes).QueryAttribute<string>("meshname");
       select_by_value      << attribute<int>   ("cycle", File.GetCycles()[timeState]);
       select_by_value      << attribute<int>   ("timelevel", 0);
       // loop over all components of the current mesh
       int nblocks = 0;
       int ngroups = 0;
       list<domain> ldomains;
       for (CarpetN5::Iterator iComps = File.Begin(select_by_identifier, select_by_value); !iComps.Done(); ++iComps, ++nblocks)
       {
          //cout << "comp: " << iComps.current_objectname() << endl;
          ldomains.push_back(*iComps);
          const int rl = ldomains.back().reflevel();
          if (ngroups < rl+1) ngroups = rl+1;
       }
       // set maximum number of reflevels
       nlevels = ngroups;
       
       // Here's the call that tells the meta-data object that we have a mesh:
       string meshname = (*iMeshes).QueryAttribute<string>("meshname").value();
       avtMeshType meshType;
       if (meshname == "Carpet-AMR") meshType = AVT_AMR_MESH;
       if (meshname == "Curvilinear") meshType = AVT_CURVILINEAR_MESH;
       
       int maxndims = 3;
       
       AddMeshToMetaData(md, meshname, meshType, NULL, nblocks, 0, maxndims, maxndims);
          
       avtMeshMetaData *mesh = (avtMeshMetaData *) md->GetMesh(meshname);
       mesh->xUnits = mesh->yUnits = mesh->zUnits = "M";
       mesh->hasSpatialExtents = false;
       mesh->blockTitle = "comps";
       mesh->blockPieceName = "comp";
       mesh->numGroups = ngroups;
       mesh->groupTitle = "reflevels";
       mesh->groupPieceName = "reflevel";
       vector<int> groupIds(nblocks);
       vector<string> blockPieceNames(nblocks);
       vector<domain> domains(nblocks);
       int i = 0;
       for (list<domain>::const_iterator ild=ldomains.begin(); ild != ldomains.end(); ++ild, ++i)
       {
          stringstream Name;
          domains[i] = *ild;
          int comp = domains[i].comp();
          int level = domains[i].reflevel();
          groupIds[i] = level; 
          Name << "ref-level=" << level << ", component=" << comp;
          blockPieceNames[i] = Name.str();
       }
       mesh->blockNames = blockPieceNames;
       md->AddGroupInformation(mesh->numGroups, nblocks, groupIds);

       if (meshname == "Carpet-AMR") domainsAMR = domains;
       if (meshname == "Curvilinear") domainsCurvilinear = domains;
       
       // loop over all variable types of the current mesh
       select_by_identifier = Nirvana::metadata(); // reset
       select_by_identifier << attribute<int>("cycle", 0);
       select_by_identifier << attribute<int>("timelevel", 0);
       select_by_identifier << attribute<string>("meshname", "");
       select_by_identifier << attribute<string>("vargrouptype", "");
       for (CarpetN5::Iterator iVarTypes = File.Begin(select_by_identifier, select_by_value); !iVarTypes.Done(); ++iVarTypes)
       {
          //cout << "vartype: " << iVarTypes.current_objectname() << endl;
          //cout << (*iVarTypes).QueryAttribute<string>("vargrouptype").value() << endl;
          
          // vector variables are defined by their groupnames
          if ((*iVarTypes).QueryAttribute<string>("vargrouptype").value() == "vector")
          {
             Nirvana::metadata select_by_identifier_v = select_by_identifier;
             Nirvana::metadata select_by_value_v = select_by_value;
             select_by_identifier_v << attribute<string>("vargroupname", "");
             select_by_value_v << (*iVarTypes).QueryAttribute<string>("vargrouptype");
             for (CarpetN5::Iterator iVars = File.Begin(select_by_identifier_v, select_by_value_v); !iVars.Done(); ++iVars) {
                //cout << "vargroup: " << iVars.current_objectname() << endl;
                // add vector variable
                int dim = 3;
                string varname = (*iVars).QueryAttribute<string>("vargroupname").value();
                string propername = varname.replace(varname.find("::"), 2, "--")+" on "+meshname;
                avtVectorMetaData *vmd = new avtVectorMetaData(propername.c_str(), meshname, AVT_NODECENT, dim);
                vmd->centering = AVT_NODECENT;
                md->Add(vmd);
                // also create variable entry
                Nirvana::metadata var_md;
                var_md << (*iVars).QueryAttribute<string>("vargrouptype");
                var_md << (*iVars).QueryAttribute<string>("vargroupname");
                vector<string> components;
                // loop over components and get the names
                Nirvana::metadata select_by_identifier_c = select_by_identifier_v;
                Nirvana::metadata select_by_value_c = select_by_value_v;
                select_by_identifier_c << attribute<string>("varname", "");
                select_by_value_c << (*iVars).QueryAttribute<string>("vargroupname");
                for (CarpetN5::Iterator iComps = File.Begin(select_by_identifier_c, select_by_value_c); !iComps.Done(); ++iComps) {
                   //cout << "component varname: " << iComps.current_objectname() << endl;
                   components.push_back((*iComps).QueryAttribute<string>("varname").value());
                }
                vectorvars[propername] = variable(var_md, components);
             }
          }
          else
          {
             // all other variables are treated as scalars and defined by their varnames
             Nirvana::metadata select_by_identifier_v = select_by_identifier;
             Nirvana::metadata select_by_value_v = select_by_value;
             select_by_identifier_v << attribute<string>("vargroupname", "");
             select_by_identifier_v << attribute<string>("varname", "");
             select_by_value_v << (*iVarTypes).QueryAttribute<string>("vargrouptype");
             for (CarpetN5::Iterator iVars = File.Begin(select_by_identifier_v, select_by_value_v); !iVars.Done(); ++iVars) {
                //cout << "varname: " << iVars.current_objectname() << endl;
                // add scalar variable
                string varname = (*iVars).QueryAttribute<string>("varname").value();
                string propername = varname.replace(varname.find("::"), 2, "--")+" on "+meshname;
                avtScalarMetaData *smd = new avtScalarMetaData(propername.c_str(), meshname, AVT_NODECENT);
                smd->centering = AVT_NODECENT;
                md->Add(smd);
                Nirvana::metadata var_md;
                // also create variable entry
                var_md << (*iVars).QueryAttribute<string>("vargrouptype");
                var_md << (*iVars).QueryAttribute<string>("vargroupname");
                var_md << (*iVars).QueryAttribute<string>("varname");
                scalarvars[propername] = variable(var_md, vector<string>(0));
             }
          }
       }
       
    }
    
    
}


// ****************************************************************************
//  Method: avtCarpetHDF5::BuildDomainAuxiliaryInfo
//
//  Purpose: Build the two data structures needed to support nesting and
//  abutting of AMR subgrids.
//  This is restricted to the Cartesian AMR grid and does not apply to the 
//  curvilinear multiblock grid.
//
//  Note: These are *never* explicitly served up to VisIt like a GetMesh or
//  GetVar call would do. Instead, we essentially 'publish' them to VisIt
//  by sticking the structures we create here into the database cache. VisIt
//  will try to look for them there when it needs them.
//
//  Programmer: Christian Reisswig
//  Creation:   Sun Feb 28 14:31:17 PDT 2010
//
// ****************************************************************************

void
avtCarpetN5FileFormat::BuildDomainAuxiliaryInfo(int timeState)
{
#ifdef MDSERVER
    return;
#endif

    vector<int> dims(3,1);

    const int num_dims = 3;
    const int num_levels = nlevels;
    const int num_patches = domainsAMR.size();

    // first, look to see if we don't already have it cached
    void_ref_ptr vrTmp = cache->GetVoidRef("any_mesh",
                                   AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
                                  timeState, -1);
    if (*vrTmp == NULL)
    {
        //
        // build the avtDomainNesting object
        //
        avtStructuredDomainNesting *dn =
            new avtStructuredDomainNesting(num_patches, num_levels);

        dn->SetNumDimensions(num_dims);

        //
        // Set refinement level ratio information
        // We hardcode this to 2:1 ratio for each direction!
        //
        vector<int> ratios(3,1);
        dn->SetLevelRefinementRatios(0, ratios);
        for (int i = 1; i < num_levels; i++)
        {
           ratios[0] = 2;
           ratios[1] = 2;
           ratios[2] = 2;
           dn->SetLevelRefinementRatios(i, ratios);
        }
        
        //
        // set each domain's level, children and logical extents
        //
        for (int i = 0; i < num_patches; i++)
        {
            vector<int> childPatches;
            float x0 = domainsAMR[i].xmin();
            float x1 = domainsAMR[i].xmax();
            float y0 = domainsAMR[i].ymin();
            float y1 = domainsAMR[i].ymax();
            float z0 = domainsAMR[i].zmin();
            float z1 = domainsAMR[i].zmax();
            for (int j = 0; j < num_patches; j++)
            {
                if (domainsAMR[j].reflevel() != domainsAMR[i].reflevel()+1)
                    continue;

                float a0 = domainsAMR[j].xmin();
                float a1 = domainsAMR[j].xmax();
                float b0 = domainsAMR[j].ymin();
                float b1 = domainsAMR[j].ymax();
                float c0 = domainsAMR[j].zmin();
                float c1 = domainsAMR[j].zmax();

                if (a0 >= x1 || x0 >= a1 ||
                    b0 >= y1 || y0 >= b1 ||
                    ((num_dims == 3) && (c0 >= z1 || z0 >= c1)))
                    continue;

                childPatches.push_back(j);
            }

            dims = domainsAMR[i].npoints();

            // fill logical extends. Note that these are CELL-indices!
            vector<int> logExts(6);
            logExts[0] = domainsAMR[i].iorigin()[0];
            logExts[1] = domainsAMR[i].iorigin()[1];
            logExts[2] = domainsAMR[i].iorigin()[2];
            logExts[3] = logExts[0] + dims[0] -2;
            logExts[4] = logExts[1] + dims[1] -2;
            //logExts[5] = num_dims == 3 ? logExts[2] + dims[2] -2 /*+ int(!isOutermostZ)*/ : 0;
            logExts[5] = logExts[2] + dims[2] -2 >= 0 ? logExts[2] + dims[2] -2 /*+ int(!isOutermostZ)*/ : 0;

            //debug1 << "logExts[0] = " << logExts[0] << ", logExts[1] = " << logExts[1] << ", logExt[2] = " << logExts[2] << endl;
            //debug1 << "logExts[3] = " << logExts[3] << ", logExts[4] = " << logExts[4] << ", logExt[5] = " << logExts[5] << endl;

            dn->SetNestingForDomain(i, domainsAMR[i].reflevel() /*-1*/,
                childPatches, logExts);
        }

        void_ref_ptr vr = void_ref_ptr(dn, avtStructuredDomainNesting::Destruct);

        cache->CacheVoidRef("any_mesh", AUXILIARY_DATA_DOMAIN_NESTING_INFORMATION,
            timeState, -1, vr);
    }
}




// ****************************************************************************
//  Method: avtCarpetN5FileFormat::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: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

vtkDataSet *
avtCarpetN5FileFormat::GetMesh(int timestate, int domain, const char *meshname)
{
   string mn(meshname);

   if (mn == "Carpet-AMR")
   {
      const int ndims = domainsAMR[domain].dim();

      vtkFloatArray* coords[3] = {0, 0, 0};
      
      // Feed AMR domain abutting to VisIt's cache
      BuildDomainAuxiliaryInfo(timestate);
      
      for (int dim = 0; dim < domainsAMR[domain].dim(); ++dim) {
        coords[dim] = vtkFloatArray::New();
        if (dim < ndims) {
          coords[dim]->SetNumberOfTuples(domainsAMR[domain].npoints()[dim]);
          float *array = (float *)coords[dim]->GetVoidPointer(0);
          for (int j=0; j < domainsAMR[domain].npoints()[dim]; j++) {
            array[j] = domainsAMR[domain].origin()[dim] +  j * domainsAMR[domain].delta()[dim];
          }
        } else {
          coords[dim]->SetNumberOfTuples(1);
          coords[dim]->SetComponent(0, 0, 0.);
        }
      }
      
      //
      // Create the vtkRectilinearGrid object and set its dimensions
      // and coordinates.
      //
      vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
      rgrid->SetDimensions(&domainsAMR[domain].npoints()[0]);
      rgrid->SetXCoordinates(coords[0]);
      coords[0]->Delete();
      rgrid->SetYCoordinates(coords[1]);
      coords[1]->Delete();
      rgrid->SetZCoordinates(coords[2]);
      coords[2]->Delete();
    
      return rgrid;
   }
   
   if (mn == "Curvilinear") 
   {   
      /*const int ndims = domainsCurvilinear[domain].dim();
      
      // This is the array that contains the coordinates
      vector<vector<double> > xyzarray(ndims);
      
      // get metadata for domain
      Nirvana::metadata md = domainCurvilinear[domain].md();
      // get metadata for coordinates
      
      // set metadata for hyperslabbing
      vector<int> offset(ndims, 0);
      for (int d=0; d < ndims; ++d)
         offset[d] = domainCurvilinear[domain].ghosts()[2*d];
      md << attribute<vector<int> >("hyperslab-offset", offset);
      md << attribute<vector<int> >("hyperslab-count", domainCurvilinear[domain].npoints());
      md << attribute<vector<int> >("hyperslab-stride", vector<int>(ndims, 1));

      vector<int> dims = domainsCurvilinear[domain].npoints();
      const int nnodes = dims[0]*dims[1] * (dims.size() > 2 ? dims[2] : 1);

      // Read the coordinates from file
      for (int i=0; i < ndims; ++i)
      {
         Nirvana::metadata md_comp = md;
         //md_comp << attribute<string>("varname", domainCurvilinear[domain].components()[i]);
         //File.GetField(md_comp, xyzarray[i]);
         assert(ntuples == xyzarray[i].size());
      } 

      //
      // Create the vtkStructuredGrid and vtkPoints objects.
      //
      vtkStructuredGrid *sgrid  = vtkStructuredGrid::New();
      vtkPoints         *points = vtkPoints::New();
      sgrid->SetPoints(points);
      sgrid->SetDimensions(&dims[0]);
      points->Delete();
      points->SetNumberOfPoints(nnodes);
      //
      // Copy the coordinate values into the vtkPoints object.
      //
      float *pts = (float *) points->GetVoidPointer(0);
      float *xc = &xyzarray[0].front();
      float *yc = &xyzarray[1].front();
      float *zc = xyzarray.size() > 2 ? &xyzarray[2].front() : NULL;
      if(ndims == 3)
      {
         for(int k = 0; k < dims[2]; ++k)
         for(int j = 0; j < dims[1]; ++j)
         for(int i = 0; i < dims[0]; ++i)
         {
            *pts++ = *xc++;
            *pts++ = *yc++;
            *pts++ = *zc++;
         }
      }
      else if(ndims == 2)
      {
         for(int j = 0; j < dims[1]; ++j)
         for(int i = 0; i < dims[0]; ++i)
         {
            *pts++ = *xc++;
            *pts++ = *yc++;
            *pts++ = 0.;
         }
      }
      
      return sgrid;*/
   }
 
   return NULL;
}


// ****************************************************************************
//  Method: avtCarpetN5FileFormat::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.
//      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: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

vtkDataArray *
avtCarpetN5FileFormat::GetVar(int timestate, int domainID, const char *varname)
{
   string str = varname;
   const bool isAMR = (str.find("Carpet-AMR") != string::npos);
   domain dom = isAMR ? domainsAMR[domainID] : domainsCurvilinear[domainID];
   
   // get metadata for domain
   Nirvana::metadata md = dom.md();
   // get metadata for variable
   md << scalarvars[str].md();
   // set metadata for hyperslabbing
   vector<int> offset(dom.dim(), 0);
   for (int d=0; d < dom.dim(); ++d)
      offset[d] = dom.ghosts()[2*d];
   md << attribute<vector<int> >("hyperslab-offset", offset);
   md << attribute<vector<int> >("hyperslab-count", dom.npoints());
   md << attribute<vector<int> >("hyperslab-stride", vector<int>(dom.dim(), 1));
   
   vector<double> data;
   File.GetField(md, data);

   int ntuples = dom.npoints()[0]
               * dom.npoints()[1]
               * (dom.npoints().size() > 2 ? dom.npoints()[2] : 1);
   
   assert(ntuples == data.size());
   
   vtkFloatArray *rv = vtkFloatArray::New();
   rv->SetNumberOfTuples(ntuples);
   for (int i = 0 ; i < ntuples ; i++) {
      rv->SetTuple1(i, data[i]);
   }
   
   return rv;
}


// ****************************************************************************
//  Method: avtCarpetN5FileFormat::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.
//      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: reisswig -- generated by xml2avt
//  Creation:   Fri Aug 27 06:08:45 PDT 2010
//
// ****************************************************************************

vtkDataArray *
avtCarpetN5FileFormat::GetVectorVar(int timestate, int domainID,const char *varname)
{

   string str = varname;
   const bool isAMR = (str.find("Carpet-AMR") != string::npos);
   domain dom = isAMR ? domainsAMR[domainID] : domainsCurvilinear[domainID];
   
   // get metadata for domain
   Nirvana::metadata md = dom.md();
   // get metadata for variable
   md << vectorvars[str].md();
   // set metadata for hyperslabbing
   vector<int> offset(dom.dim(), 0);
   for (int d=0; d < dom.dim(); ++d)
      offset[d] = dom.ghosts()[2*d];
   md << attribute<vector<int> >("hyperslab-offset", offset);
   md << attribute<vector<int> >("hyperslab-count", dom.npoints());
   md << attribute<vector<int> >("hyperslab-stride", vector<int>(dom.dim(), 1));
   
   const int ncomps = vectorvars[str].dim();  // This is the rank of the vector - typically 2 or 3.
   
   int ntuples = dom.npoints()[0]
               * dom.npoints()[1]
               * dom.npoints()[2];
   
   // read in each component separately
   vector<vector<double> > data(ncomps);
   for (int i=0; i < ncomps; ++i)
   {
      Nirvana::metadata md_comp = md;
      md_comp << attribute<string>("varname", vectorvars[str].components()[i]);
      File.GetField(md_comp, data[i]);
      assert(ntuples == data[i].size());
   }
   
   vtkFloatArray *rv = vtkFloatArray::New();
   const int ucomps = (ncomps == 2 ? 3 : ncomps);
   rv->SetNumberOfComponents(ucomps);
   rv->SetNumberOfTuples(ntuples);
   float *one_entry = new float[ucomps];
   for (int i = 0 ; i < ntuples ; i++)
   {
      int j;
      for (j = 0 ; j < ncomps ; j++)
           one_entry[j] = data[j][i];
      for (j = ncomps ; j < ucomps ; j++)
           one_entry[j] = 0.;
      rv->SetTuple(i, one_entry); 
   }
   
   delete [] one_entry;
   return rv;
}