Commit 9fa8f225 authored by George Zagaris's avatar George Zagaris
Browse files

ENH:Move internal enzo reader to a separate header

This commit moves the internal enzo reader to a
separate header. The reason for doing this is
primarily to allow the re-use of the internal enzo
reader for the corresponding particles reader.
parent 0572d68b
......@@ -43,1360 +43,9 @@
#include <vtkstd/string>
#include <cassert>
//==============================================================================
#include "vtkAMREnzoReaderInternal.hpp"
/*****************************************************************************
*
* Copyright (c) 2000 - 2009, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-400124
* All rights reserved.
*
* This file was adapted from the VisIt Enzo reader (avtEnzoFileFormat). 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.
*
*****************************************************************************/
#define ENZO_READER_SLASH_CHAR '\\'
#define ENZO_READER_SLASH_STRING "\\"
const int ENZO_READER_BUFFER_SIZE = 4096;
static char ENZO_READER_STRING[ ENZO_READER_BUFFER_SIZE ];
// ----------------------------------------------------------------------------
// Functions for Parsing File Names
// ----------------------------------------------------------------------------
static const char * GetEnzoMajorFileName( const char * path, int & start )
{
start = 0;
if (path == 0)
{
strcpy( ENZO_READER_STRING, "." );
return ENZO_READER_STRING;
}
else
if ( *path == '\0' )
{
strcpy( ENZO_READER_STRING, "." );
return ENZO_READER_STRING;
}
else
{
// find end of path string
int n = 0;
while ( ( path[n] != '\0' ) && ( n < ENZO_READER_BUFFER_SIZE ) )
{
n ++;
}
// deal with string too large
if ( n == ENZO_READER_BUFFER_SIZE )
{
strcpy( ENZO_READER_STRING, "." );
return ENZO_READER_STRING;
}
// backup, skipping over all trailing ENZO_READER_SLASH_CHAR chars
int j = n-1;
while ( ( j >= 0 ) && ( path[j] == ENZO_READER_SLASH_CHAR ) )
{
j --;
}
// deal with string consisting of all ENZO_READER_SLASH_CHAR chars
if ( j == -1 )
{
start = -1;
strcpy( ENZO_READER_STRING, ENZO_READER_SLASH_STRING );
return ENZO_READER_STRING;
}
// backup to just after next ENZO_READER_SLASH_CHAR char
int i = j-1;
while ( ( i >= 0 ) && ( path[i] != ENZO_READER_SLASH_CHAR ) )
{
i --;
}
i ++;
start = i;
// build the return string
int k;
for ( k = 0; k < j - i + 1; k ++ )
{
ENZO_READER_STRING[k] = path[ i + k ];
}
ENZO_READER_STRING[k] = '\0';
return ENZO_READER_STRING;
}
}
const char * GetEnzoMajorFileName( const char * path )
{
// int dummy1;
// return GetEnzoMajorFileName( path, dummy1 );
std::vector< std::string > vpath;
vtksys::SystemTools::SplitPath( path,vpath );
assert( vpath.size() >= 1);
return( vpath[ vpath.size()-1 ].c_str() );
}
const char * GetEnzoDirectory( const char * path )
{
int start;
GetEnzoMajorFileName( path, start );
std::string mydir = vtksys::SystemTools::GetFilenamePath( std::string(path) );
return mydir.c_str( );
}
// ----------------------------------------------------------------------------
// Class vtkEnzoReaderBlock (begin)
// ----------------------------------------------------------------------------
class vtkEnzoReaderBlock
{
public:
vtkEnzoReaderBlock() { this->Init(); }
~vtkEnzoReaderBlock() { this->Init(); }
int Index;
int Level;
int ParentId;
vtkstd::vector< int > ChildrenIds;
int MinParentWiseIds[3];
int MaxParentWiseIds[3];
int MinLevelBasedIds[3];
int MaxLevelBasedIds[3];
int NumberOfParticles;
int NumberOfDimensions;
int BlockCellDimensions[3];
int BlockNodeDimensions[3];
double MinBounds[3];
double MaxBounds[3];
double SubdivisionRatio[3];
vtkstd::string BlockFileName;
vtkstd::string ParticleFileName;
void Init()
{
this->BlockFileName = "";
this->ParticleFileName = "";
this->Index = -1;
this->Level = -1;
this->ParentId = -1;
this->ChildrenIds.clear();
this->NumberOfParticles = 0;
this->NumberOfDimensions = 0;
this->MinParentWiseIds[0] =
this->MinParentWiseIds[1] =
this->MinParentWiseIds[2] =
this->MaxParentWiseIds[0] =
this->MaxParentWiseIds[1] =
this->MaxParentWiseIds[2] = -1;
this->MinLevelBasedIds[0] =
this->MinLevelBasedIds[1] =
this->MinLevelBasedIds[2] =
this->MaxLevelBasedIds[0] =
this->MaxLevelBasedIds[1] =
this->MaxLevelBasedIds[2] = -1;
this->BlockCellDimensions[0] =
this->BlockCellDimensions[1] =
this->BlockCellDimensions[2] =
this->BlockNodeDimensions[0] =
this->BlockNodeDimensions[1] =
this->BlockNodeDimensions[2] = 0;
this->MinBounds[0] =
this->MinBounds[1] =
this->MinBounds[2] = VTK_DOUBLE_MAX;
this->MaxBounds[0] =
this->MaxBounds[1] =
this->MaxBounds[2] =-VTK_DOUBLE_MAX;
this->SubdivisionRatio[0] =
this->SubdivisionRatio[1] =
this->SubdivisionRatio[2] = 1.0;
}
void GetParentWiseIds( vtkstd::vector< vtkEnzoReaderBlock > & blocks );
void GetLevelBasedIds( vtkstd::vector< vtkEnzoReaderBlock > & blocks );
};
// ----------------------------------------------------------------------------
// get the bounding (cell) Ids of this block in terms of its parent block's
// sub-division resolution (indexing is limited to the scope of the parent)
void vtkEnzoReaderBlock::GetParentWiseIds
( vtkstd::vector< vtkEnzoReaderBlock > & blocks )
{
if ( this->ParentId != 0 )
{
// the parent is not the root and then we need to determine the offset
// (in terms of the number of parent divisions / cells) of the current
// block's beginning / ending position relative to the parent block's
// beginning position
vtkEnzoReaderBlock & parent = blocks[ this->ParentId ];
this->MinParentWiseIds[0] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[0]
* ( this->MinBounds[0] - parent.MinBounds[0] )
/ ( parent.MaxBounds[0] - parent.MinBounds[0] ) );
this->MaxParentWiseIds[0] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[0]
* ( this->MaxBounds[0] - parent.MinBounds[0] )
/ ( parent.MaxBounds[0] - parent.MinBounds[0] ) );
this->MinParentWiseIds[1] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[1]
* ( this->MinBounds[1] - parent.MinBounds[1] )
/ ( parent.MaxBounds[1] - parent.MinBounds[1] ) );
this->MaxParentWiseIds[1] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[1]
* ( this->MaxBounds[1] - parent.MinBounds[1] )
/ ( parent.MaxBounds[1] - parent.MinBounds[1] ) );
if ( this->NumberOfDimensions == 3 )
{
this->MinParentWiseIds[2] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[2]
* ( this->MinBounds[2] - parent.MinBounds[2] )
/ ( parent.MaxBounds[2] - parent.MinBounds[2] ) );
this->MaxParentWiseIds[2] = static_cast < int >
( 0.5 + parent.BlockCellDimensions[2]
* ( this->MaxBounds[2] - parent.MinBounds[2] )
/ ( parent.MaxBounds[2] - parent.MinBounds[2] ) );
}
else
{
this->MinParentWiseIds[2] = 0;
this->MaxParentWiseIds[2] = 0;
}
// the ratio for mapping two parent-wise Ids to 0 and
// this->BlockCellDimension[i],
// respectively, while the same region is covered
this->SubdivisionRatio[0] = static_cast < double >
( this->BlockCellDimensions[0] ) /
( this->MaxParentWiseIds[0] - this->MinParentWiseIds[0] );
this->SubdivisionRatio[1] = static_cast < double >
( this->BlockCellDimensions[1] ) /
( this->MaxParentWiseIds[1] - this->MinParentWiseIds[1] );
if ( this->NumberOfDimensions == 3 )
{
this->SubdivisionRatio[2] = static_cast < double >
( this->BlockCellDimensions[2] ) /
( this->MaxParentWiseIds[2] - this->MinParentWiseIds[2] );
}
else
{
this->SubdivisionRatio[2] = 1.0;
}
}
else
{
// Now that the parent is the root, it can not provide cell-dimensions
// information (BlockCellDimensions[0 .. 2]) directly, as the above does.
// Fortunately we can obtain it according to the spatial ratio of the
// child block (the current one) to the parent (root) and the child block's
// cell-dimensions information. This derivation is based on the definition
// of 'level' that all children blocks at the same level (e.g., the current
// block and its siblings) share the same sub-division ratio relative to
// their parent (the root herein).
vtkEnzoReaderBlock & block0 = blocks[0];
double xRatio = ( this->MaxBounds[0] - this->MinBounds[0] ) /
( block0.MaxBounds[0] - block0.MinBounds[0] );
this->MinParentWiseIds[0] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[0] / xRatio ) // parent's dim
* ( this->MinBounds[0] - block0.MinBounds[0] )
/ ( block0.MaxBounds[0] - block0.MinBounds[0] )
);
this->MaxParentWiseIds[0] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[0] / xRatio )
* ( this->MaxBounds[0] - block0.MinBounds[0] )
/ ( block0.MaxBounds[0] - block0.MinBounds[0] )
);
double yRatio = ( this->MaxBounds[1] - this->MinBounds[1] ) /
( block0.MaxBounds[1] - block0.MinBounds[1] );
this->MinParentWiseIds[1] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[1] / yRatio )
* ( this->MinBounds[1] - block0.MinBounds[1] )
/ ( block0.MaxBounds[1] - block0.MinBounds[1] )
);
this->MaxParentWiseIds[1] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[1] / yRatio )
* ( this->MaxBounds[1] - block0.MinBounds[1] )
/ ( block0.MaxBounds[1] - block0.MinBounds[1] )
);
if ( this->NumberOfDimensions == 3 )
{
double zRatio = ( this->MaxBounds[2] - this->MinBounds[2] ) /
( block0.MaxBounds[2] - block0.MinBounds[2] );
this->MinParentWiseIds[2] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[2] / zRatio )
* ( this->MinBounds[2] - block0.MinBounds[2] )
/ ( block0.MaxBounds[2] - block0.MinBounds[2] )
);
this->MaxParentWiseIds[2] = static_cast < int >
( 0.5 + ( this->BlockCellDimensions[2] / zRatio )
* ( this->MaxBounds[2] - block0.MinBounds[2] )
/ ( block0.MaxBounds[2] - block0.MinBounds[2] )
);
}
else
{
this->MinParentWiseIds[2] = 0;
this->MaxParentWiseIds[2] = 0;
}
this->SubdivisionRatio[0] = 1.0;
this->SubdivisionRatio[1] = 1.0;
this->SubdivisionRatio[2] = 1.0;
}
}
// ----------------------------------------------------------------------------
// determine the bounding (cell) Ids of this block in terms of the sub-division
// resolution of the level at which its parent block lies (indexing includes
// all siblings of this block --- those sibling blocks beyond the scope of the
// parent of this block)
void vtkEnzoReaderBlock::GetLevelBasedIds
( vtkstd::vector< vtkEnzoReaderBlock > & blocks )
{
// note that this function is invoked from the root in a top-down manner
// and the parent-wise Ids have been determined in advance
if ( this->ParentId != 0 )
{
// the parent is not the root and therefore we need to exploit the level-
// based Ids of the parent, of which the shifted verson is multiplied by
// the refinement ratio
vtkEnzoReaderBlock & parent = blocks[ this->ParentId ];
this->MinLevelBasedIds[0] = static_cast < int >
( ( parent.MinLevelBasedIds[0]
+ this->MinParentWiseIds[0]
) * this->SubdivisionRatio[0]
);
this->MinLevelBasedIds[1] = static_cast < int >
( ( parent.MinLevelBasedIds[1]
+ this->MinParentWiseIds[1]
) * this->SubdivisionRatio[1]
);
this->MinLevelBasedIds[2] = static_cast < int >
( ( parent.MinLevelBasedIds[2]
+ this->MinParentWiseIds[2]
) * this->SubdivisionRatio[2]
);
this->MaxLevelBasedIds[0] = static_cast < int >
( ( parent.MinLevelBasedIds[0]
+ this->MaxParentWiseIds[0]
) * this->SubdivisionRatio[0]
);
this->MaxLevelBasedIds[1] = static_cast < int >
( ( parent.MinLevelBasedIds[1]
+ this->MaxParentWiseIds[1]
) * this->SubdivisionRatio[1]
);
this->MaxLevelBasedIds[2] = static_cast < int >
( ( parent.MinLevelBasedIds[2]
+ this->MaxParentWiseIds[2]
) * this->SubdivisionRatio[2]
);
}
else
{
// now that the parent is the root, the parent-wise Ids
// are just the level-based Ids
this->MinLevelBasedIds[0] = this->MinParentWiseIds[0];
this->MinLevelBasedIds[1] = this->MinParentWiseIds[1];
this->MinLevelBasedIds[2] = this->MinParentWiseIds[2];
this->MaxLevelBasedIds[0] = this->MaxParentWiseIds[0];
this->MaxLevelBasedIds[1] = this->MaxParentWiseIds[1];
this->MaxLevelBasedIds[2] = this->MaxParentWiseIds[2];
}
}
// ----------------------------------------------------------------------------
// Class vtkEnzoReaderBlock ( end )
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
// Class vtkEnzoReaderInternal (begin)
// ----------------------------------------------------------------------------
class vtkEnzoReaderInternal
{
public:
vtkEnzoReaderInternal( vtkAMREnzoReader * reader )
{ this->Init(); this->TheReader = reader; }
vtkEnzoReaderInternal() { this->Init(); }
~vtkEnzoReaderInternal() { this->ReleaseDataArray();
this->Init();
}
// number of all vtkDataSet (vtkImageData / vtkRectilinearGrid / vtkPolyData)
// objects that have been SUCCESSFULLY extracted and inserted to the output
// vtkMultiBlockDataSet (including rectilinear blocks and particle sets)
int NumberOfMultiBlocks;
int NumberOfDimensions;
int NumberOfLevels;
int NumberOfBlocks;
int ReferenceBlock;
int CycleIndex;
char * FileName;
double DataTime;
vtkDataArray * DataArray;
vtkAMREnzoReader * TheReader;
vtkstd::string DirectoryName;
vtkstd::string MajorFileName;
vtkstd::string BoundaryFileName;
vtkstd::string HierarchyFileName;
vtkstd::vector< vtkstd::string > BlockAttributeNames;
vtkstd::vector< vtkstd::string > ParticleAttributeNames;
vtkstd::vector< vtkstd::string > TracerParticleAttributeNames;
vtkstd::vector< vtkEnzoReaderBlock > Blocks;
void Init()
{
this->DataTime = 0.0;
this->FileName = NULL;
this->TheReader = NULL;
this->DataArray = NULL;
this->CycleIndex = 0;
this->ReferenceBlock = 0;
this->NumberOfBlocks = 0;
this->NumberOfLevels = 0;
this->NumberOfDimensions = 0;
this->NumberOfMultiBlocks = 0;
this->DirectoryName = "";
this->MajorFileName = "";
this->BoundaryFileName = "";
this->HierarchyFileName = "";
this->Blocks.clear();
this->BlockAttributeNames.clear();
this->ParticleAttributeNames.clear();
this->TracerParticleAttributeNames.clear();
}
void ReleaseDataArray()
{
if ( this->DataArray )
{
this->DataArray->Delete();
this->DataArray = NULL;
}
}
void SetFileName( char * fileName ) { this->FileName = fileName; }
void ReadMetaData();
void GetAttributeNames();
void CheckAttributeNames();
void ReadBlockStructures();
void ReadGeneralParameters();
void DetermineRootBoundingBox();
int LoadAttribute( const char *attribute, int blockIdx );
int GetBlockAttribute(
const char* attribute, int blockIdx, vtkDataSet* pDataSet );
};
int vtkEnzoReaderInternal::GetBlockAttribute(
const char *atribute, int blockIdx, vtkDataSet *pDataSet )
{
// this function must be called by GetBlock( ... )
this->ReadMetaData();
if ( atribute == NULL || blockIdx < 0 ||
pDataSet == NULL || blockIdx >= this->NumberOfBlocks )
{
return 0;
}
// try obtaining the attribute and attaching it to the grid as a cell data
// NOTE: the 'equal' comparison below is MUST because in some cases (such
// as cosmological datasets) not all rectilinear blocks contain an assumably
// common block attribute. This is the case where such a block attribute
// may result from a particles-to-cells interpolation process. Thus those
// blocks with particles (e.g., the reference one with the fewest cells)
// contain it as a block attribute, whereas others without particles just
// do not contain this block attribute.
int succeded = 0;
if ( this->LoadAttribute( atribute, blockIdx ) &&
( pDataSet->GetNumberOfCells() ==
this->DataArray->GetNumberOfTuples() )
)
{
succeded = 1;
pDataSet->GetCellData()->AddArray( this->DataArray );
this->ReleaseDataArray();
}
return succeded;
}
int vtkEnzoReaderInternal::LoadAttribute( const char *atribute, int blockIdx )
{
// TODO: implement this
// called by GetBlockAttribute( ... ) or GetParticlesAttribute()
this->ReadMetaData();
if ( atribute == NULL || blockIdx < 0 ||
blockIdx >= this->NumberOfBlocks )
{
return 0;
}
// this->Internal->Blocks includes a pseudo block --- the root as block #0
blockIdx ++;
// currently only the HDF5 file format is supported
vtkstd::string blckFile = this->Blocks[ blockIdx ].BlockFileName;
hid_t fileIndx = H5Fopen( blckFile.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if( fileIndx < 0 )
{
return 0;
}
// retrieve the contents of the root directory to look for a group
// corresponding to the target block, if available, open that group
int blckIndx;
char blckName[65];
int objIndex;
hsize_t numbObjs;
hid_t rootIndx = H5Gopen( fileIndx, "/" );
H5Gget_num_objs( rootIndx, &numbObjs );
for ( objIndex=0; objIndex < static_cast < int >( numbObjs ); objIndex ++ )
{
int type = H5Gget_objtype_by_idx( rootIndx, objIndex );
if ( type == H5G_GROUP )
{
H5Gget_objname_by_idx( rootIndx, objIndex, blckName, 64 );
if ( ( sscanf( blckName, "Grid%d", &blckIndx ) == 1 ) &&
( (blckIndx == blockIdx) || (blckIndx == blockIdx+1) ) )
{
// located the target block
rootIndx = H5Gopen( rootIndx, blckName );
break;
}
}
} // END for all objects
// disable error messages while looking for the attribute (if any) name
// and enable error messages when it is done
void * pContext = NULL;
H5E_auto_t erorFunc;
H5Eget_auto( &erorFunc, &pContext );
H5Eset_auto( NULL, NULL );
hid_t attrIndx = H5Dopen( rootIndx, atribute );
H5Eset_auto( erorFunc, pContext );
pContext = NULL;
// check if the data attribute exists
if ( attrIndx < 0 )
{
vtkGenericWarningMacro(
"Attribute (" << atribute << ") data does not exist in file "
<< blckFile.c_str() );
H5Gclose( rootIndx );
H5Fclose( fileIndx );
return 0;
}
// get cell dimensions and the valid number
hsize_t cellDims[3];
hid_t spaceIdx = H5Dget_space( attrIndx );
H5Sget_simple_extent_dims( spaceIdx, cellDims, NULL );
hsize_t numbDims = H5Sget_simple_extent_ndims( spaceIdx );
// number of attribute tuples = number of cells (or particles)
int numTupls = 0;
switch( numbDims )
{
case 1: