/*****************************************************************************
 *
 * Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
 * Produced at the Lawrence Livermore National Laboratory
 * All rights reserved.
 *
 * This file is part of VisIt. For details, see http://www.llnl.gov/visit/. 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 materials provided with the distribution.
 *  - Neither the name of the UC/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  THE  REGENTS  OF  THE  UNIVERSITY OF
 * CALIFORNIA, 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.
 *
 *****************************************************************************/

// ************************************************************************* //
//                            avtlataFileFormat.C                           //
// ************************************************************************* //

#include <avtlataFileFormat.h>

#include <LmlReader.h>
#include <LataJournal.h>

#include <avtDatabaseMetaData.h>
#include <avtGhostData.h>
#include <DebugStream.h>
#include <Expression.h>
#include <InvalidVariableException.h>

#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkFloatArray.h>
#include <vtkInformation.h>
#include <vtkIntArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkStructuredGrid.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnstructuredGrid.h>

#include <string>
#include <fstream>
#include <iostream>

// ****************************************************************************
//  Method: avtlata constructor
//
//  Programmer: fauchet -- generated by xml2avt
//
// ****************************************************************************

avtlataFileFormat::avtlataFileFormat(const char *filename)
  : avtMTMDFileFormat(filename)
{
  debug1 << "avtlataFileFormat constructor " << filename << endl;
  try {
    set_Journal_level(0);

    LataOptions opt;
    LataOptions::extract_path_basename(filename, opt.path_prefix, opt.basename);
    opt.dual_mesh = true;
    opt.faces_mesh = true;
    opt.regularize = 2;
    opt.regularize_tolerance = 1e-7;
    opt.user_fields_=true;
    read_any_format_options(filename, opt);
    debug1 << "avtlataFileFormat: initializing filter" << endl;
    // Read the source file to the lata database
    read_any_format(filename, opt.path_prefix, lata_db_);
    filter_.initialize(opt, lata_db_);
  }
  catch (LataDBError err) {
    cerr << "Error in LataFilter::initialize " << filename << " " << err.describe() << endl;
    throw;
  }
}

avtlataFileFormat::~avtlataFileFormat()
{
}

// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: fauchet -- generated by xml2avt
//
// ****************************************************************************

int
avtlataFileFormat::GetNTimesteps(void)
{
  int n;
  try {
    n = filter_.get_nb_timesteps();
    // Timestep 0 contains global definitions.
    // If we have "real" timesteps, do not show timestep 0
    if (n > 1)
      n--;
  }
  catch (LataDBError err) {
    cerr << "Error in getntimesteps " << filename << " " << err.describe() << endl;
    throw;
  }
  return n;
}

void avtlataFileFormat::GetTimes(std::vector<double>& times) 
{ 
  int n;
  try {
    n = filter_.get_nb_timesteps();
    if (n == 1)
      times.push_back(0.);
    else
      for (int i = 1; i < n; i++)
        times.push_back(filter_.get_timestep(i));
  }
  catch (LataDBError err) {
    cerr << "Error in gettimes " << filename << " " << err.describe() << endl;
    throw;
  }
  return;
}

// ****************************************************************************
//  Method: avtlataFileFormat::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: fauchet -- generated by xml2avt
//
// ****************************************************************************

void
avtlataFileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avtlataFileFormat::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: fauchet -- generated by xml2avt
//
// ****************************************************************************

void
avtlataFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md, int timeState)
{
  try {
  debug1 << "avtlataFileFormat::PopulateDatabaseMetaData : " 
         << filename << " " << timeState << endl;

  const char *suffix_vector_names[] = { "_X", "_Y", "_Z" };
  const char *suffix_vector_expr[] = { "[0]", "[1]", "[2]" };

  const Noms geoms = filter_.get_exportable_geometry_names();

  for (int i_geom = 0; i_geom < geoms.size(); i_geom++) {
    debug1 << " Domain : " << geoms[i_geom] << endl;
    const LataGeometryMetaData data = filter_.get_geometry_metadata(geoms[i_geom]);

    avtMeshType mt = AVT_UNSTRUCTURED_MESH;
    if (data.internal_name_.finit_par("IJK"))
      {
        mt = AVT_RECTILINEAR_MESH;
      }
    int block_origin = 0;
    int topo_dim;
    switch(data.element_type_) {
    case Domain::point:     topo_dim = 0; mt = AVT_POINT_MESH; break;
    case Domain::line:      topo_dim = 1; break;
    case Domain::triangle: 
    case Domain::quadri:    topo_dim = 2; break;
    case Domain::tetra:
    case Domain::prism6:
    case Domain::polyedre:
    case Domain::hexa:      topo_dim = 3; break;
    default:
      topo_dim = 3; ///TODO: this should be an error in default case!
      EXCEPTION1(InvalidVariableException, 
                 "avtlataFileFormat::PopulateDatabaseMetaData error: unknown element type");
      cerr << "avtlataFileFormat::PopulateDatabaseMetaData error: unknown element type" << endl;
//      exit(-1);
    }

    int mesh_faces=0;
    if  (data.internal_name_.finit_par("_centerfaces"))
      {
        //cerr<<"la "<<data.internal_name_<<endl;
        mesh_faces=1;
      }
    double *extents = NULL;
    const std::string geom_name(data.displayed_name_);
    AddMeshToMetaData(md, geom_name, mt, extents, data.nblocks_, block_origin,
                      data.dimension_, topo_dim);
    mesh_username_.add(data.displayed_name_);
    mesh_latafilter_name_.add(data.internal_name_);

    Field_UNames fields = filter_.get_exportable_field_unames(geoms[i_geom]);

    for (int i_field = 0; i_field < fields.size(); i_field++) {
      const LataFieldMetaData data2 = filter_.get_field_metadata(fields[i_field]);
      avtCentering cent;
      switch (data2.localisation_) {
      case LataField_base::ELEM: cent = AVT_ZONECENT; break;
      case LataField_base::SOM:  cent = AVT_NODECENT; break;
      default:
        // Do not export fields that cannot be shown
        continue;
      }

      // Take localisation of source field
      Nom loc = data2.source_localisation_;
      std::string varname(data2.name_);
      varname += "_";
      varname += loc;
      varname += "_";
      varname += geom_name;
      if (data2.nb_components_ == 1) {
        // Scalar field
        // We append the geometry name to the component name:
        register_fieldname(varname.c_str(), fields[i_field], 0);
        AddScalarVarToMetaData(md, varname, geom_name, cent);
      } else if (data2.is_vector_ && data2.nb_components_ == data.dimension_) {
        // Vector field
        register_fieldname(varname.c_str(), fields[i_field], -1);
        AddVectorVarToMetaData(md, varname, geom_name, cent, data2.nb_components_);
        if (mesh_faces==0)
          {
        std::string n;
        for (entier i = 0; i < data2.nb_components_; i++) {
          Expression v;
          n = data2.name_;
          n += suffix_vector_names[i];
          n += "_";
          n += loc;
          n += "_";
          n += geom_name;
          v.SetName(n);
          n = varname;
          n += suffix_vector_expr[i];
          v.SetDefinition(n);
          v.SetType(Expression::ScalarMeshVar);
          md->AddExpression(&v);
        }
        Expression norme_v;
        n = "norme_";
        n += varname;
        norme_v.SetName(n);
        n = "magnitude(";
        n += varname;
        n += ")";
        norme_v.SetDefinition(n);
        norme_v.SetType(Expression::ScalarMeshVar);
        md->AddExpression(&norme_v);
          }
      } else {
        // Multiscalar field
        // I chose to postfix the varname with the component name, perhaps not the best choice.
        if (mesh_faces==0)
          {
        for (entier i_compo = 0; i_compo < data2.nb_components_; i_compo++) {
          std::string varname2(data2.name_);
          varname2 += "_";
          if (data2.component_names_.size() == data2.nb_components_) {
            varname2 += data2.component_names_[i_compo];
          } else {
            Nom n(i_compo);
            varname2 += n;
          }
          varname2 += "_";
          varname2 += loc;
          varname2 += "_";
          varname2 += geom_name;
          register_fieldname(varname2.c_str(), fields[i_field], i_compo);
          AddScalarVarToMetaData(md, varname2, geom_name, cent);
        }
          }
      }
    }
  }
  debug1 << "End avtlataFileFormat::PopulateDatabaseMetaData" << endl;
  }
  catch (LataDBError err) {
    cerr << "Error in PopulateDatabaseMetaData " << err.describe() << endl;
    throw;
  }
}

void
avtlataFileFormat::register_fieldname(const char *visit_name, const Field_UName & uname, int component)
{
  if (field_username_.rang(visit_name) >= 0) {
    cerr << "Error in avtlataFileFormat::register_fieldname: duplicate field name " << visit_name << endl;
    cerr << "Ignoring field" << endl;
    return;
  }
  field_username_.add(visit_name);
  field_uname_.add(uname);
  field_component_.add(component);
}

void
avtlataFileFormat::register_meshname(const char *visit_name, const char *latafilter_name)
{
  if (mesh_username_.rang(visit_name) >= 0) {
    cerr << "Error in avtlataFileFormat::register_meshname: duplicate name " << visit_name << endl;
    cerr << "Ignoring mesh" << endl;
    return;    
  }
  mesh_username_.add(visit_name);
  mesh_latafilter_name_.add(latafilter_name);
}

// ****************************************************************************
//  Method: avtlataFileFormat::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: fauchet -- generated by xml2avt
//
// ****************************************************************************

vtkDataSet *
avtlataFileFormat::GetMesh(int timestate, int block, const char *meshname)
{
  vtkDataSet *return_value = 0;
  try {
  debug1 << " avtlataFileFormat::GetMesh ts=" << timestate 
         << " block=" << block 
         << " meshname=" << meshname << endl;

  // We have real timesteps in the database, add one to timestep index:
  if (filter_.get_nb_timesteps() > 1)
    timestate++;

  const entier index = mesh_username_.rang(meshname);
  if (index < 0) {
    cerr << "internal error in avtlataFileFormat::GetMesh: name " << meshname << " not found" << endl;
//    exit(-1);
  }
  Domain_Id id(mesh_latafilter_name_[index], timestate, block);
  const Domain & geometry = filter_.get_geometry(id);

  const DomainUnstructured * geom_ptr = dynamic_cast<const DomainUnstructured*>(&geometry);
  const DomainIJK          * ijk_ptr = dynamic_cast<const DomainIJK*>(&geometry);

  if (geom_ptr) {
    const DomainUnstructured & geom = *geom_ptr;

    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    vtkPoints *points = vtkPoints::New();
    const FloatTab & pos = geom.nodes_;
    const int nnodes = pos.dimension(0);
    const int dim3 = pos.dimension(1) == 3;
    points->SetNumberOfPoints(nnodes);
    float* pts = (float *) points->GetVoidPointer(0);
    int jl=0;
    int i;
    for (i = 0; i < nnodes; i++) {
      pts[jl]   = pos(i,0);
      pts[jl+1] = pos(i,1);
      pts[jl+2] = dim3 ? pos(i,2) : 0.;
      jl+=3;
    }
    ugrid->SetPoints(points);
    points->Delete();
    
    const IntTab & conn = geom.elements_;
    const int ncells = conn.dimension(0);
    int nverts = conn.dimension(1);
    
    int type_cell;
    switch (geom.elt_type_) {
    case Domain::point:
      type_cell=VTK_VERTEX; 
      if (ncells == 0) 
        nverts = 1;
      break;
    case Domain::line:
      type_cell=VTK_LINE;
      break;
    case Domain::triangle:
      type_cell=VTK_TRIANGLE;
      break;
    case Domain::quadri:
      type_cell=VTK_QUAD;
      break;
    case Domain::tetra:
      type_cell=VTK_TETRA;
      break;
    case Domain::prism6:
      type_cell=VTK_WEDGE;
      break;
    case Domain::hexa:
      type_cell=VTK_HEXAHEDRON;
      break;
    case Domain::polyedre:
      type_cell=VTK_CONVEX_POINT_SET;
      break;
    default:
      type_cell=-1;
      cerr<<"avtlataFileFormat::GetMesh unknown elt type "<<endl;
//      exit(-1);
      break;
    }
    vtkIdType *verts = new vtkIdType[nverts];
    if (type_cell == VTK_VERTEX && ncells == 0) {
      // Cells are implicit. Create them:
      ugrid->Allocate(nnodes);
      for (i = 0; i < nnodes; i++) {
        verts[0] = i;
        ugrid->InsertNextCell(type_cell, nverts, verts);
      }
    } else {
      ugrid->Allocate(ncells);
      for (i = 0; i < ncells; i++) {
        if (type_cell==VTK_QUAD) {
          // Nodes order is different in visit than in trio_u
          verts[0]=conn(i,0);
          verts[1]=conn(i,1);
          verts[2]=conn(i,3);
          verts[3]=conn(i,2);
        } else if (type_cell==VTK_HEXAHEDRON) {
          // Nodes order is different in visit than in trio_u
          verts[0]=conn(i,0);
          verts[1]=conn(i,1);
          verts[2]=conn(i,3);
          verts[3]=conn(i,2);
          verts[4]=conn(i,4);
          verts[5]=conn(i,5);
          verts[6]=conn(i,7);
          verts[7]=conn(i,6);
        } else if (type_cell==VTK_CONVEX_POINT_SET) {
          int nverts_loc=nverts;
          for (int j = 0; j < nverts; j++) 
            {
              verts[j] = conn(i,j);
             
              if (verts[j]==-1)
                {  
                  nverts_loc=j; 
                  break;
                }    
            }
          int  nb_som_max_to_regularize=0;
          if (filter_.get_options().regularize_polyedre!=0)
            {
              nb_som_max_to_regularize=8;
              if (filter_.get_options().regularize_polyedre==-1)
                  nb_som_max_to_regularize=32000;
            }
          if ((nb_som_max_to_regularize>=6) && (nverts_loc==6))
            ugrid->InsertNextCell(VTK_WEDGE, nverts_loc, verts);
          else if ((nb_som_max_to_regularize>=12)&&(nverts_loc==12))
            ugrid->InsertNextCell(VTK_HEXAGONAL_PRISM, nverts_loc, verts);
          else if ((nb_som_max_to_regularize>=8)&&(nverts_loc==8))
            {
                        // Nodes order is different in visit than in trio_u
              verts[0]=conn(i,0);
              verts[1]=conn(i,1);
              verts[2]=conn(i,3);
              verts[3]=conn(i,2);
              verts[4]=conn(i,4);
              verts[5]=conn(i,5);
              verts[6]=conn(i,7);
              verts[7]=conn(i,6);
              ugrid->InsertNextCell(VTK_HEXAHEDRON, nverts_loc, verts);
               
            }
          else
            ugrid->InsertNextCell(type_cell, nverts_loc, verts);
        }
        else {
          for (int j = 0; j < nverts; j++) 
            verts[j] = conn(i,j);
        }
        if (type_cell!=VTK_CONVEX_POINT_SET)          
         
          ugrid->InsertNextCell(type_cell, nverts, verts);
      }
    }
    delete [] verts;
    verts = 0;
    // Declare ghost elements:
    const int n = geom.nb_virt_items(LataField_base::ELEM);
    if (n > 0) {
      unsigned char realVal = 0;
      unsigned char ghost   = 0; // Sera modifie par AddGhostZoneType
      avtGhostData::AddGhostZoneType(ghost, DUPLICATED_ZONE_INTERNAL_TO_PROBLEM);
      vtkUnsignedCharArray *ghostcells = vtkUnsignedCharArray::New();
      ghostcells->SetName("avtGhostZones");
      ghostcells->SetNumberOfTuples(ncells);
      unsigned char *dat = (unsigned char *) ghostcells->GetVoidPointer(0);
      for (i = 0; i < ncells - n; i++)
        dat[i] = realVal;
      for (i = ncells - n; i < ncells; i++)
        dat[i] = ghost;
      ugrid->GetCellData()->AddArray(ghostcells);
      ugrid->GetInformation()->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), 0);
      ghostcells->Delete();
    }
    return_value = ugrid;
    
  } else if (ijk_ptr) {
    const DomainIJK & geom = *ijk_ptr;

    // Maillage regulier : on transmet la grille ijk
    vtkRectilinearGrid *sgrid = vtkRectilinearGrid::New();

    const int dim = geom.coord_.size();
    ArrOfInt ncoord(3, 1);
    int i;
    for (i = 0; i < dim; i++) 
      ncoord[i] = geom.coord_[i].size_array();
    sgrid->SetDimensions(ncoord[0], ncoord[1], ncoord[2]);

    for (i = 0; i < 3; i++) {
      float *data;
      vtkFloatArray *c;
      c = vtkFloatArray::New();
      const int n = ncoord[i];
      c->SetNumberOfTuples(n);
      data = (float *) c->GetVoidPointer(0);
      if (i < dim) {
        const ArrOfFloat & coord = geom.coord_[i];
        for (int j = 0; j < n; j++)
          data[j] = coord[j];
      } else {
        data[0] = 0.;
      }
      switch(i) {
      case 0: sgrid->SetXCoordinates(c); break;
      case 1: sgrid->SetYCoordinates(c); break;
      case 2: sgrid->SetZCoordinates(c); break;
      default: ;
      }
      c->Delete();
    }
    // Create "invalid cells" data (GettingDataIntoVisit.pdf, page 136)
    // and "ghost cells"
    const int n = geom.invalid_connections_.size_array();
    if (n > 0 || geom.virtual_layer_begin_ || geom.virtual_layer_end_) {
      const int ncells = geom.nb_elements();
      unsigned char realVal = 0;
      unsigned char invalid = 0; // Sera modifie par AddGhostZoneType
      unsigned char ghost   = 0;
      avtGhostData::AddGhostZoneType(invalid, ZONE_NOT_APPLICABLE_TO_PROBLEM);
      avtGhostData::AddGhostZoneType(ghost, DUPLICATED_ZONE_INTERNAL_TO_PROBLEM);
      vtkUnsignedCharArray *ghostcells = vtkUnsignedCharArray::New();
      ghostcells->SetName("avtGhostZones");
      ghostcells->SetNumberOfTuples(ncells);
      unsigned char *dat = (unsigned char *) ghostcells->GetVoidPointer(0);

      for (i = 0; i < ncells; i++)
          dat[i] = realVal;

      if (n > 0) {
        // invalid cells
        for (i = 0; i < ncells; i++) {
          if (geom.invalid_connections_[i])
            dat[i] = invalid;
        }
      }
      
      // ghost cells
      entier ij = 1;
      for (i = 0; i < dim-1; i++)
        ij *= ncoord[i]-1;
      if (geom.virtual_layer_begin_) {
        // first layer of cells is ghost
        for (i = 0; i < ij * geom.virtual_layer_begin_; i++)
          dat[i] += ghost;
      }
      if (geom.virtual_layer_end_) {
        // last layer of cells is ghost
        for (i = ncells - ij * geom.virtual_layer_end_; i < ncells; i++)
          dat[i] += ghost;
      }
      
      sgrid->GetCellData()->AddArray(ghostcells);
      sgrid->GetInformation()->Set(
        vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), 0);
      ghostcells->Delete();
    }

    return_value = sgrid;
  } else {
    cerr << "Error in avtlataFileFormat::GetMesh: unknown geometry type" << endl;
//    exit(-1);
  }

  filter_.release_geometry(geometry);
  }
  catch (LataDBError err) {
    cerr << "Error in getmesh " << timestate << " " << block << " " << meshname << " " << err.describe() << endl;
    throw;
  }

  return return_value;
}


// ****************************************************************************
//  Method: avtlataFileFormat::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: fauchet -- generated by xml2avt
//
// ****************************************************************************

vtkDataArray *
avtlataFileFormat::GetVar(int timestate, int block, const char *varname)
{
  vtkDataArray * return_value = 0;
  try {
  debug1 << "Getvar time:" << timestate
         << " block:" << block
         << " varname:" << varname << endl;

  if (filter_.get_nb_timesteps() > 1)
    timestate++;

  Field_UName field_uname;
  int component;
  get_field_info_from_visitname(varname, field_uname, component);

  if (component < 0) {
    cerr << "Error: avtlataFileFormat::GetVar called for vector field" << endl;
//    exit(-1);
  }

  Field_Id id(field_uname, timestate, block);

  const LataField_base & field = filter_.get_field(id);

  const Field<FloatTab> * float_field_ptr = dynamic_cast<const Field<FloatTab>*>(&field);
  const Field<IntTab> * int_field_ptr = dynamic_cast<const Field<IntTab>*>(&field);

  if (float_field_ptr) {
    vtkFloatArray *rv = vtkFloatArray::New();
    const Field<FloatTab> & fld = *float_field_ptr;
    const FloatTab & values = fld.data_;
    int ntuples = values.dimension(0);
    rv->SetNumberOfTuples(ntuples);
    float * data = rv->GetPointer(0);
    for (int i = 0; i < ntuples; i++) 
      data[i] = values(i, component);
    return_value = rv;
  } else if (int_field_ptr) {
    vtkIntArray *rv = vtkIntArray::New();
    const Field<IntTab> & fld = *int_field_ptr;
    const IntTab & values = fld.data_;
    int ntuples = values.dimension(0);
    rv->SetNumberOfTuples(ntuples);
    int * data = rv->GetPointer(0);
    for (int i = 0; i < ntuples; i++) 
      data[i] = values(i, component);    
    return_value = rv;
  } else {
    cerr << "Error in avtlataFileFormat::GetVar: unknown data type" << endl;
//    exit(-1);
  }
  filter_.release_field(field);
  }
  catch (LataDBError err) {
    cerr << "Error in getvar " << timestate << " " << block << " " << varname << " " << err.describe() << endl;
    throw;
  }
  return return_value; 
}


// ****************************************************************************
//  Method: avtlataFileFormat::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: fauchet -- generated by xml2avt
//
// ****************************************************************************

vtkDataArray *
avtlataFileFormat::GetVectorVar(int timestate, int block, const char *varname)
{
  vtkDataArray * return_value = 0;
  try {
  debug1 << "Getvectorvar time:" << timestate
         << " block:" << block
         << " varname:" << varname << endl;

  if (filter_.get_nb_timesteps() > 1)
    timestate++;

  Field_UName field_uname;
  int component;
  get_field_info_from_visitname(varname, field_uname, component);

  if (component >= 0) {
    cerr << "Error: avtlataFileFormat::GetVectorVar called for scalar field" << endl;
//    exit(-1);
  }

  Field_Id id(field_uname, timestate, block);

  const LataField_base & field = filter_.get_field(id);

  const Field<FloatTab> * float_field_ptr = dynamic_cast<const Field<FloatTab>*>(&field);
  const Field<IntTab> * int_field_ptr = dynamic_cast<const Field<IntTab>*>(&field);

  if (float_field_ptr) {
    vtkFloatArray *rv = vtkFloatArray::New();
    const Field<FloatTab> & fld = *float_field_ptr;
    const FloatTab & values = fld.data_;
    int ntuples = values.dimension(0);
    int dim = values.dimension(1);
    rv->SetNumberOfComponents(3);
    rv->SetNumberOfTuples(ntuples);
    float* data= rv->WritePointer(0,3*ntuples);
    for (int i = 0; i < ntuples; i++)
      for (int j = 0; j < 3; j++)
        data[i*3+j] = (j<dim) ? values(i, j) : 0.;
    return_value = rv;
  } else if (int_field_ptr) {
    vtkIntArray *rv = vtkIntArray::New();
    const Field<IntTab> & fld = *int_field_ptr;
    const IntTab & values = fld.data_;
    int ntuples = values.dimension(0);
    int dim = values.dimension(1);
    rv->SetNumberOfComponents(3);
    rv->SetNumberOfTuples(ntuples);
    int* data= rv->WritePointer(0,3*ntuples);
    for (int i = 0; i < ntuples; i++) 
      for (int j = 0; j < 3; j++)
        data[i*3+j] = (j<dim) ? values(i, j) : 0;
    return_value = rv;
  } else {
    cerr << "Error in avtlataFileFormat::GetVectorVar: unknown data type" << endl;
//    exit(-1);
  }
  filter_.release_field(field);
  }
  catch (LataDBError err) {
    cerr << "Error in getvectorvar " << timestate << " " << block << " " << varname << " " << err.describe() << endl;
    throw;
  }
  return return_value; 
}

void avtlataFileFormat::get_field_info_from_visitname(const char *varname, Field_UName & uname, int & component) const
{
  const int k = field_username_.rang(varname);
  if (k < 0) {
    cerr << "Error in avtlataFileFormat::get_field_info_from_visitname: field " << varname << " not found" << endl;
//    exit(-1);
  }
  uname = field_uname_[k];
  component = field_component_[k];
}
