Commits

Douglas Rudd  committed 4b6c73a

Updated version of artio, replaced _artio_reader, other modifications to artio frontend

  • Participants
  • Parent commits f5c46f2

Comments (0)

Files changed (23)

File yt/frontends/artio/_artio_caller.pyx

 from libc.stdlib cimport malloc, free
 import  data_structures  
 
-ROOT_LEVEL = 1 #snlroot this has to be 1 currently 
-
-cdef struct artio_context_struct: 
-    int comm
-ctypedef artio_context_struct * artio_context
-cdef artio_context artio_context_global
-
-cdef extern from "sfc.h":
-    int sfc_index( int coords[3], int num_root_grid_refinements ) 
-    int sfc_index0( int ix, int iy, int iz, int num_root_grid_refinements ) 
-
 cdef extern from "artio.h":
-    ctypedef struct artio_file_struct "artio_file_struct" :
+    ctypedef struct artio_fileset_handle "artio_fileset" :
         pass
-    ctypedef artio_file_struct *artio_file "artio_file"
-    ctypedef struct artio_grid_file_struct "artio_grid_file_struct" :
+    ctypedef struct artio_context :
         pass
-    ctypedef artio_grid_file_struct *artio_grid_file "artio_grid_file"
-
-
-    ctypedef struct artio_context_struct "artio_context_struct" :
-        pass
-    ctypedef artio_context_struct *artio_context "artio_context"
-
-
-    cdef artio_context artio_context_global
-    
+    cdef extern artio_context *artio_context_global 
 
     # open modes
     cdef int ARTIO_OPEN_HEADER "ARTIO_OPEN_HEADER"
     cdef int ARTIO_READ_REFINED "ARTIO_READ_REFINED"
     cdef int ARTIO_READ_ALL "ARTIO_READ_ALL"
     cdef int ARTIO_READ_REFINED_NOT_ROOT "ARTIO_READ_REFINED_NOT_ROOT"
-    
+    cdef int ARTIO_RETURN_CELLS "ARTIO_RETURN_CELLS"
+    cdef int ARTIO_RETURN_OCTS "ARTIO_RETURN_OCTS"
 
     # errors
     cdef int ARTIO_SUCCESS "ARTIO_SUCCESS"
     cdef int ARTIO_ERR_MEMORY_ALLOCATION "ARTIO_ERR_MEMORY_ALLOCATION"
 
-    artio_file artio_fileset_open(char *file_prefix, int type, artio_context context )
-    int artio_fileset_close( artio_file handle )
-    int artio_fileset_open_particle( artio_file handle )
-    int artio_fileset_open_grid(artio_file handle) 
-    int artio_fileset_close_grid(artio_file handle) 
+    artio_fileset_handle *artio_fileset_open(char *file_prefix, int type, artio_context *context )
+    int artio_fileset_close( artio_fileset_handle *handle )
+    int artio_fileset_open_particle( artio_fileset_handle *handle )
+    int artio_fileset_open_grid(artio_fileset_handle *handle) 
+    int artio_fileset_close_grid(artio_fileset_handle *handle) 
 
     # parameter functions
-    int artio_parameter_iterate( artio_file handle, char *key, int *type, int *length )
-    int artio_parameter_get_int_array(artio_file handle, char * key, int length, int32_t *values)
-    int artio_parameter_get_float_array(artio_file handle, char * key, int length, float *values)
-    int artio_parameter_get_long_array(artio_file handle, char * key, int length, int64_t *values)
-    int artio_parameter_get_double_array(artio_file handle, char * key, int length, double *values)
-    int artio_parameter_get_string_array(artio_file handle, char * key, int length, char **values, int max_length)
-    int artio_grid_read_root_cell_end(artio_file handle)
-    int artio_grid_read_root_nocts(artio_file handle, int64_t sfc,\
-                                  float *variables, int32_t *num_oct_levels, int32_t *num_octs_per_level)
-    int artio_grid_cache_sfc_range(artio_file handle, int64_t start, int64_t end)
+    int artio_parameter_iterate( artio_fileset_handle *handle, char *key, int *type, int *length )
+    int artio_parameter_get_int_array(artio_fileset_handle *handle, char * key, int length, int32_t *values)
+    int artio_parameter_get_float_array(artio_fileset_handle *handle, char * key, int length, float *values)
+    int artio_parameter_get_long_array(artio_fileset_handle *handle, char * key, int length, int64_t *values)
+    int artio_parameter_get_double_array(artio_fileset_handle *handle, char * key, int length, double *values)
+    int artio_parameter_get_string_array(artio_fileset_handle *handle, char * key, int length, char **values, int max_length)
 
-    ctypedef void (* GridCallBackYTPos)(float * variables, int level, int refined, int64_t isfc, double pos[3], void *)
-    int artio_grid_read_sfc_range_ytpos(artio_file handle,\
-                int64_t sfc_min, int64_t sfc_max,\
-                int min_level_to_read, int max_level_to_read,\
-                int options,\
-                GridCallBackYTPos callback, void *pyobject)
+    int artio_grid_cache_sfc_range(artio_fileset_handle *handle, int64_t start, int64_t end)
+    int artio_grid_clear_sfc_cache( artio_fileset_handle *handle ) 
 
-    ctypedef void (* GridCallBackYT)(float * variables, int level, int refined, int64_t isfc, void *)
-    int artio_grid_read_sfc_range_yt(artio_file handle,\
-                int64_t sfc_min, int64_t sfc_max,\
-                int min_level_to_read, int max_level_to_read,\
-                int options,\
-                GridCallBackYT callback, void *pyobject)
+    int artio_grid_read_root_cell_begin(artio_fileset_handle *handle, int64_t sfc, 
+        double *pos, float *variables,
+        int *num_tree_levels, int *num_octs_per_level)
+    int artio_grid_read_root_cell_end(artio_fileset_handle *handle)
 
-    ctypedef void (* GridCallBack)(float * variables, int level, int refined,int64_t isfc)
-    int artio_grid_read_sfc_range(artio_file handle,\
-                int64_t sfc_min, int64_t sfc_max,\
-                int min_level_to_read, int max_level_to_read,\
-                int options,\
-                GridCallBack callback)
-    
-                                 
-    ctypedef void (* ParticleCallBack)(int pid, double *primary_variables,\
-                                           double *secondary_variables, int species,\
-                                           int subspecies, int isfc, void *pyobject )
-    int artio_particle_read_sfc_range_yt(artio_file handle,\
-                                          int64_t sfc_min, int64_t sfc_max,\
-                                          int64_t species_min, int64_t species_max,\
-                                          ParticleCallBack callback, void *pyobject)
-########## wrappers calling c
-cdef extern from "artio.c":
-    artio_file artio_fileset_open(char * file_prefix, int type, artio_context context) 
+    int artio_grid_read_level_begin(artio_fileset_handle *handle, int level )
+    int artio_grid_read_level_end(artio_fileset_handle *handle)
 
+    int artio_grid_read_oct(artio_fileset_handle *handle, double *pos, 
+            float *variables, int *refined)
 
-cdef class read_parameters_artio : 
+    int artio_grid_count_octs_in_sfc_range(artio_fileset_handle *handle,
+            int64_t start, int64_t end, int64_t *num_octs)
+
+def check_artio_status(status, fname="[unknown]"):
+    if status!=ARTIO_SUCCESS :
+        callername = sys._getframe().f_code.co_name
+        nline = sys._getframe().f_lineno
+        print 'failure with status', status, 'in function',fname,'from caller', callername, nline 
+        sys.exit(1)
+
+cdef class artio_fileset :
     cdef public object parameters 
-    cdef artio_file handle
+    cdef object oct_handler
+    cdef artio_fileset_handle *handle
+    cdef int64_t num_root_cells
+    cdef int64_t sfc_min, sfc_max
+    cdef public int num_grid
 
-    def __init__(self, char *file_prefix, int artio_type) :
+    # grid attributes
+    cdef int min_level, max_level
+    cdef int num_grid_variables
+
+    cdef int cpu
+    cdef int domain_id
+
+    def __init__(self, char *file_prefix) :
+        cdef int artio_type = ARTIO_OPEN_HEADER
+        cdef int64_t num_root
+
         self.handle = artio_fileset_open( file_prefix, artio_type, artio_context_global ) 
         self.read_parameters()
-        artio_fileset_close(self.handle)  
+        print 'print parameters in caller.pyx',self.parameters
+        print 'done reading header parameters'
+
+        self.num_root_cells = self.parameters['num_root_cells'][0]
+        self.num_grid = 1
+        num_root = self.num_root_cells
+        while num_root > 1 :
+            self.num_grid <<= 1
+            num_root >>= 3
+ 
+        # dhr - add grid detection code 
+        status = artio_fileset_open_grid( self.handle )
+        check_artio_status(status)
+  
+
+        # grid stuff
+        self.oct_handler=None
+
+        self.min_level = 0
+        self.max_level = self.parameters['grid_max_level'][0]
+
+        #snl FIX: the sfc values used should come from "subset" and describe the domain for chunking
+        # note the root level method may force chunking to be done on 0-level ytocts 
+        self.sfc_min = 0
+        self.sfc_max = self.parameters['grid_file_sfc_index'][1]-1
+        self.num_grid_variables = self.parameters['num_grid_variables'][0]
+
+        # these should be fixed to something meaningful
+        self.cpu = 0
+        self.domain_id = 0
 
     def read_parameters(self) :
         cdef char key[64]
 
             self.parameters[key] = parameter
 
-def check_artio_status(status, fname):
-    callername = sys._getframe().f_code.co_name
-    nline = sys._getframe().f_lineno
-    if status!=ARTIO_SUCCESS :
-        print 'failure with status', status, 'in function',fname,'from caller', callername, nline 
-        sys.exit(1)
-cdef class artio_fileset :
-    cdef public object parameters 
-    cdef artio_file handle
-    def __init__(self, char *file_prefix) :
-        cdef int artio_type = ARTIO_OPEN_HEADER
-        self.handle = artio_fileset_open( file_prefix, artio_type, artio_context_global ) 
-        d = read_parameters_artio(file_prefix, artio_type)
-        self.parameters = {}
-        self.parameters = d.parameters
-        print 'print parameters in caller.pyx',self.parameters
-        print 'done reading header parameters'
+    def count_refined_octs(self) :
+        cdef int64_t num_total_octs = 0
 
-cdef class artio_fileset_grid :
-    cdef public object parameters
-    cdef artio_file handle
-    cdef public file_prefix
+        # this only works if this domain includes all cells!
+        if self.parameters.has_key("num_octs_per_level") :
+            return self.parameters["num_octs_per_level"].sum()
 
-    def __init__(self, char *file_prefix) :
-        cdef int artio_type = ARTIO_OPEN_GRID
-        self.handle = artio_fileset_open( file_prefix, artio_type, artio_context_global ) 
-        artio_fileset_open_grid( self.handle ) 
-        d = read_parameters_artio(file_prefix, artio_type)
-        self.parameters = {}
-        self.parameters = d.parameters
-        self.file_prefix = file_prefix
-        print 'done reading grid parameters'
+        status = artio_grid_count_octs_in_sfc_range( self.handle, 
+                self.sfc_min, self.sfc_max, &num_total_octs )
+        check_artio_status(status) 
 
-#snl: subclass some of this?
-class artio_grid_routines(object) : 
-    def __init__(self, param_handle) :
-        self.oct_handler=None
-        self.source = None
-        self.art_level_count = None
-
-        self.min_level_to_read = 0
-        self.max_level_to_read = param_handle.parameters['grid_max_level'][0]
-        #snl FIX: the sfc values used should come from "subset" and describe the domain for chunking
-        # note the root level method may force chunking to be done on 0-level ytocts 
-        self.sfc_min = 0   
-        self.sfc_max = param_handle.parameters['grid_file_sfc_index'][1]-1
-        self.num_grid_variables = param_handle.parameters['num_grid_variables'][0]
-        self.num_root_grid_refinements = int( np.log10( 
-                param_handle.parameters["num_root_cells"][0]+0.5) / 
-                                              (3*np.log10(2)))
-        self.param_handle=param_handle
-        self.ng=1 
-        self.cpu=0
-        self.domain_id=0
-
-        self.grid_variable_labels=param_handle.parameters['grid_variable_labels']
+        # add octs for root cells
+        num_total_octs += (self.sfc_max-self.sfc_min+1)/8
  
-        # the dictionary from artio variable names to yt 
-        #snl FIX: this should in fields.py ?
-        self.arttoyt_label_var = {}
-        for label in self.grid_variable_labels :
-            if   label == 'HVAR_GAS_DENSITY' :
-                self.arttoyt_label_var[label] = 'Density'
-            elif label == 'HVAR_GAS_ENERGY' :
-                self.arttoyt_label_var[label] = 'TotalGasEnergy'
-            elif label == 'HVAR_INTERNAL_ENERGY' :
-                self.arttoyt_label_var[label] = 'GasEnergy'
-            elif label == 'HVAR_PRESSURE' :
-                self.arttoyt_label_var[label] = 'Pressure'
-            elif label == 'HVAR_MOMENTUM_X' :
-                self.arttoyt_label_var[label] = 'XMomentumDensity'
-            elif label == 'HVAR_MOMENTUM_Y' :
-                self.arttoyt_label_var[label] = 'YMomentumDensity'
-            elif label == 'HVAR_MOMENTUM_Z' :
-                self.arttoyt_label_var[label] = 'ZMomentumDensity'
-            elif label == 'HVAR_GAMMA' :
-                self.arttoyt_label_var[label] = 'Gamma'
-            elif label == 'HVAR_METAL_DENSITY_II' :
-                self.arttoyt_label_var[label] = 'MetalDensitySNIa'
-            elif label == 'HVAR_METAL_DENSITY_Ia' :
-                self.arttoyt_label_var[label] = 'MetalDensitySNII'
-            elif label == 'VAR_POTENTIAL' :
-                self.arttoyt_label_var[label] = 'Potential'
-            elif label == 'VAR_POTENTIAL_HYDRO' :
-                self.arttoyt_label_var[label] = 'PotentialHydro'
-            else :
-                self.arttoyt_label_var[label] = label
-        print 'arttoyt_label_var (in caller.pyx):', self.arttoyt_label_var
-
-        self.label_index = {}
-        self.matched_fieldnames = []
-        #snl FIX not sure about file handling:
-        # self.handle = <object> handle #<------seg faults
-        # and you dont bother to close all of these handles
-        
-    def count_refined_octs(self) :
-        cdef int min_level_to_read = self.min_level_to_read
-        cdef int max_level_to_read = self.max_level_to_read
-        cdef int64_t sfc_min = self.sfc_min
-        cdef int64_t sfc_max = self.sfc_max
-        cdef num_grid_variables = self.num_grid_variables
-        
-        cdef float * variables  
-        cdef int32_t * num_oct_levels 
-        cdef int32_t * num_octs_per_level 
-        
-        length = num_grid_variables * 8
-        variables = <float *>malloc(length*sizeof(float))
-        num_oct_levels = <int32_t *>malloc(1*sizeof(int32_t))
-        length = max_level_to_read
-        num_octs_per_level = <int32_t *>malloc(length*sizeof(int32_t))
-        
-        cdef artio_file handle
-        handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                     ARTIO_OPEN_GRID, artio_context_global ) 
-        artio_fileset_open_grid( handle ) 
-        
-        cdef int64_t num_total_octs =0
-        n_levels = max_level_to_read - min_level_to_read + 1
-        art_level_count = np.zeros(n_levels, dtype='int64')
-
-        status = artio_grid_cache_sfc_range(handle, sfc_min, sfc_max)
-        check_artio_status(status, artio_grid_routines.__name__)
-        for sfc in xrange(sfc_min,sfc_max):
-            status = artio_grid_read_root_nocts(handle, sfc,
-                                                variables, num_oct_levels,
-                                                num_octs_per_level)
-            check_artio_status(status, artio_grid_routines.__name__)
-            noct_levels = num_oct_levels[0]
-            count_level_octs = {}          
-            count_level_octs = [ num_octs_per_level[i] for i in xrange(noct_levels) ]
-            for level in xrange(noct_levels) : 
-                art_level_count[level+1] += count_level_octs[level]*8 
-            num_sfc_octs = sum(count_level_octs)
-            num_total_octs += num_sfc_octs
-            status = artio_grid_read_root_cell_end(handle)
-            check_artio_status(status, artio_grid_routines.__name__)
-
-        if ROOT_LEVEL == 1 :    
-            #add root-level octs to the count
-            art_level_count[0] = (self.sfc_max+1)
-            num_total_octs += (self.sfc_max+1)/8 
-        self.art_level_count = art_level_count 
-        print '_artio_caller.pyx num_total_octs (too big compared to what pos counts?)', num_total_octs
         return num_total_octs
 
-    def grid_pos_fill(self, oct_handler, domain_dimensions1D) :
+    def grid_pos_fill(self, oct_handler) :
         ''' adds all refined octs and a new array of ghost octs for  
         the "fully refined octs" at level=-1 in ART or 0 in yt convention 
         so that root level can consist of children
         '''
+        cdef int64_t sfc
+        cdef int level
+        cdef int num_oct_levels
+        cdef int *num_octs_per_level
+        cdef double dpos[3]
+        cdef int oct_count
+
+        pos = np.empty((1,3), dtype='float64')
+
         print 'start filling oct positions'
         self.oct_handler = oct_handler
-        self.oct_count=0
-        # fill the root grid yt-only octs 
-        if ROOT_LEVEL == 1 :    
-            pos = np.empty((1,3), dtype='float64')
-            self.domain_dimensions = domain_dimensions1D
-            for iz in  range(self.domain_dimensions/2) :
-                for iy in  range(self.domain_dimensions/2) :
-                    for ix in range(self.domain_dimensions/2) :
-                        pos[0,0]=ix*2+1
-                        pos[0,1]=iy*2+1
-                        pos[0,2]=iz*2+1
-                        level=0
-                        self.oct_count += 1
-                        self.oct_handler.add(self.cpu+1, level, 
-                                             self.ng, pos, self.domain_id)
-            ###################################
+
+        oct_count = 0
+        level = 0
+        for iz in range(self.num_grid/2) :
+            pos[0,2] = iz*2+1
+            for iy in range(self.num_grid/2) :
+                pos[0,1] = iy*2+1
+                for ix in range(self.num_grid/2) :
+                    pos[0,0]=ix*2+1
+                    self.oct_handler.add(self.cpu+1, level, self.num_grid, pos, self.domain_id) 
+                    oct_count += 1
+
         # Now do real ART octs
         print 'start filling oct positions children'
-        cdef artio_file handle
-        handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                     ARTIO_OPEN_GRID, artio_context_global ) 
-        status = artio_grid_read_sfc_range_ytpos(handle,\
-                    self.sfc_min, self.sfc_max,\
-                    self.min_level_to_read, self.max_level_to_read,\
-                    ARTIO_READ_REFINED,\
-                    wrap_oct_pos_callback, <void*>self) 
-        check_artio_status(status, artio_grid_routines.__name__)
-        artio_fileset_close(handle) 
-        print 'done filling oct positions; allocated octs:', self.oct_count
-        # snl FIX assert oct_count matches num octs elsewhere
+        status = artio_grid_cache_sfc_range( self.handle, self.sfc_min, self.sfc_max )
+        check_artio_status(status) 
+
+        num_octs_per_level = <int *>malloc(self.max_level*sizeof(int))
+
+        for sfc in range( self.sfc_min, self.sfc_max+1 ) :
+            status = artio_grid_read_root_cell_begin( self.handle, sfc, 
+                dpos, NULL, &num_oct_levels, num_octs_per_level )
+            check_artio_status(status) 
+
+            for level in range(num_oct_levels) :
+                status = artio_grid_read_level_begin( self.handle, level+1 )
+                check_artio_status(status) 
+
+                for oct in range(num_octs_per_level[level]) :
+                    status = artio_grid_read_oct( self.handle, dpos, NULL, NULL )
+                    check_artio_status(status) 
+ 
+                    pos[0,0] = dpos[0]
+                    pos[0,1] = dpos[1]
+                    pos[0,2] = dpos[2]
+                    oct_count += self.oct_handler.add(self.cpu+1, level+1, self.num_grid, pos, self.domain_id )
+
+                status = artio_grid_read_level_end( self.handle )
+                check_artio_status(status) 
+
+            status = artio_grid_read_root_cell_end( self.handle )
+            check_artio_status(status) 
+        
+        status = artio_grid_clear_sfc_cache( self.handle )
+        check_artio_status(status)
+
+        free(num_octs_per_level) 
+
+        print 'done filling oct positions', oct_count
+
     def grid_var_fill(self, source, fields):
-        print 'start filling grid vars the root grid fill takes too long...'
-        self.source = source
-        i=-1
-        for artlabel in self.grid_variable_labels :
-            label = self.arttoyt_label_var[artlabel]
-            i=i+1
-            for field in fields : 
-                if field == label :
-                    print 'match in fields?',artlabel, field,label, fields
-                    self.label_index[field]=i
-                    self.matched_fieldnames.append(field)
-        print 'matched fields:',self.matched_fieldnames
-        print 'art index of matched fields',self.label_index
+        cdef int num_oct_levels
+        cdef int *num_octs_per_level
+        cdef float *variables
+        cdef int status
+        cdef int root_oct, child, order
+        cdef int ix, iy, iz
+        cdef int cx, cy, cz
+        cdef int count
+        cdef double dpos[3]
 
+        print "Field list:", fields
+ 
+        # translate fields from ARTIO names to indices
+        if not all(f in self.parameters['grid_variable_labels'] for f in fields) :
+            print "Asked for a variable that doesn't exist!"
+            sys.exit(1)
+        field_order = dict([(f,self.parameters['grid_variable_labels'].index(f)) for f in fields])
+
+        print "field order: ", field_order
+ 
+        status = artio_grid_cache_sfc_range( self.handle, self.sfc_min, self.sfc_max )
+        check_artio_status(status) 
+
+        num_octs_per_level = <int *>malloc(self.max_level*sizeof(int))
+        variables = <float *>malloc(8*self.num_grid_variables*sizeof(float))
+
+        #art_order = {}
+        #for ix in range(2) :
+        #    for iy in range(2) :
+        #        for iz in range(2) :
+        #        #    art_order[ix+2*(iy+2*iz)] = iz+2*(iy+2*ix)
+        #            art_order[iz+2*(iy+2*ix)] = ix+2*(iy+2*iz)
+        #
+        #for i in range(8) :
+        #    print i, art_order[i]
+
+        count = self.num_root_cells
+        seen = [False for i in range(self.num_root_cells)]
+
+        for sfc in range( self.sfc_min, self.sfc_max+1 ) :
+            status = artio_grid_read_root_cell_begin( self.handle, sfc, 
+                    dpos, variables, &num_oct_levels, num_octs_per_level )
+            check_artio_status(status) 
+
+            ix = (int)(dpos[0]-0.5) / 2
+            iy = (int)(dpos[1]-0.5) / 2
+            iz = (int)(dpos[2]-0.5) / 2
+
+            cx = 0 if dpos[0] < (2*ix + 1) else 1
+            cy = 0 if dpos[1] < (2*iy + 1) else 1
+            cz = 0 if dpos[2] < (2*iz + 1) else 1
+            
+            root_oct = iz+(self.num_grid/2)*(iy+(self.num_grid/2)*ix)
+            child = cz+2*(cy+2*cx)
+            order = 8*root_oct + child
+
+            assert( root_oct < self.num_root_cells / 8 )
+            assert( child >= 0 and child < 8 )
+            assert( order >= 0 and order < self.num_root_cells )
+
+            assert( not seen[order] )
+            seen[order] = True
+
+            for f in fields :
+                source[f][order] = variables[field_order[f]]
+ 
+            for level in range(num_oct_levels) :
+                status = artio_grid_read_level_begin( self.handle, level+1 )
+                check_artio_status(status) 
+
+                for oct in range(num_octs_per_level[level]) :
+                    status = artio_grid_read_oct( self.handle, NULL, variables, NULL )
+                    check_artio_status(status) 
+
+                    for child in range(8) :
+                        for f in fields :
+                            source[f][count] = variables[self.num_grid_variables*child+field_order[f]]
+                        count += 1
+ 
+                status = artio_grid_read_level_end( self.handle )
+                check_artio_status(status) 
+
+            status = artio_grid_read_root_cell_end( self.handle )
+            check_artio_status(status) 
         
-#        i=-1
-#        for field in fields : 
-#            i+=1
-#            print field
-#            self.matched_fieldnames.append(field)
-#            self.label_index[field] = i
-#        print 'quitting in caller'
-#        sys.exit(1)
-        self.count=0
-        #fill grid positions with values
-        cdef artio_file handle
-        
-        if ROOT_LEVEL == 1 :    
-            coords = np.empty(3, dtype='int64')
-            child = np.empty(3, dtype='int64')
-            if len(self.label_index) > 0 :
-                # start with root run over root-level yt-only octs 
-                # where yt-octs live on level 0 
-                handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                             ARTIO_OPEN_GRID, artio_context_global ) 
-                for iz in  range(self.domain_dimensions/2) :
-                    for iy in  range(self.domain_dimensions/2) :
-                        for ix in range(self.domain_dimensions/2) :
-                            coords[0]=ix*2
-                            coords[1]=iy*2
-                            coords[2]=iz*2
-                            for k in range(2):
-                                for j in range(2):
-                                    for i in range(2):
-                                        child[0] = coords[0]+i
-                                        child[1] = coords[1]+j
-                                        child[2] = coords[2]+k
-                                        isfc = sfc_index0(child[0], child[1], child[2], self.num_root_grid_refinements)
-                                        status = artio_grid_read_sfc_range_yt(
-                                            handle, isfc, isfc,\
-                                                0, 0,\
-                                                ARTIO_READ_ALL,\
-                                                wrap_cell_var_callback,\
-                                                <void*>self
-                                            ) #only reads root
-                                        check_artio_status(status, artio_grid_routines.__name__)
-                artio_fileset_close(handle) 
-                print 'snlroot done with root var fill'
-            #now run through oct children
-            handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                         ARTIO_OPEN_GRID, artio_context_global ) 
-            status = artio_grid_read_sfc_range_yt(
-                handle, self.sfc_min, self.sfc_max,\
-                    self.min_level_to_read+1, self.max_level_to_read,\
-                    ARTIO_READ_ALL,\
-                    wrap_cell_var_callback,\
-                    <void*>self
-                ) #only octs!
-            check_artio_status(status, artio_grid_routines.__name__)
-            artio_fileset_close(handle) 
-        print 'done buffering variables'
-    def oct_pos_callback(self, level, refined, isfc, pos):
-#        print 'callerpos ',self.oct_count*8,pos[0,0],pos[0,1],pos[0,2],vars, level
-        self.oct_count += 1
-        self.oct_handler.add(self.cpu+1, level-self.min_level_to_read, 
-                             self.ng, pos, self.domain_id)
-    def cell_var_callback(self, level, refined, ichild, cell_var):
-        for field in self.matched_fieldnames : 
-            self.source[field][self.count] = cell_var[self.label_index[field]] 
-#        if (ichild == 0) and (level>0):
-#            print 'callervar ',self.count,cell_var[0],level
-        self.count += 1
- 
+        status = artio_grid_clear_sfc_cache( self.handle )
+        check_artio_status(status)
 
-###### callbacks (e.g. https://github.com/cython/cython/tree/master/Demos/callback) ######
-cdef void wrap_oct_pos_callback(float *variables, int level, int refined, 
-                                 int64_t isfc, double *pos, void *pyobject):
-    position = np.empty((1, 3), dtype='float64')
-    position[0,0] = pos[0]
-    position[0,1] = pos[1]
-    position[0,2] = pos[2]
-    # add one to level because in yt, octs live on the same level as their children
-    # 0-level ART octs do not exist in memory (the ART root cells are not children)
-    ytlevel = level+1 
-    artioroutines = <object>pyobject
-    #    print '_artio_caller.pyx:octpositionsandvalues',ytlevel, pos[0],pos[1],pos[2],level,variables[0]
-    artioroutines.oct_pos_callback(ytlevel, refined, isfc, position) #variables[0]
+        free(num_octs_per_level) 
+        free(variables)
 
-cdef void wrap_cell_var_callback(float *variables, int level, int refined, 
-                                 int64_t ichild, void *pyobject):
-    artioroutines = <object>pyobject
-    cell_var={}
-    cell_var = [ variables[i] for i in range(artioroutines.num_grid_variables) ]
-    #    if level > 0:
-    #        print '_artio_caller.pyx:sourcecallvalue', level, cell_var
-    artioroutines.cell_var_callback(level, refined, ichild, cell_var)
+        print 'done filling oct variables', count
 
-class artio_particle_routines(object) : 
-    def __init__(self, param_handle) :
-
-        self.min_level_to_read = 0
-        self.max_level_to_read = param_handle.parameters['grid_max_level'][0]
-        # sfc values should come from a subset object and describe the chunked domains
-        self.sfc_min = 0   
-        self.sfc_max = param_handle.parameters['grid_file_sfc_index'][1]-1
-        self.param_handle=param_handle
-        self.ng=1 #RISM
-        self.cpu=0 #RISM
-        self.domain_id=0 #RISM
-
-        # FIX: too much is hardcoded here, much of this should probably live in fields.py
-        # choose particle types (species labels), species (self.num_particle_species) 
-        # and fields (variable_labels) to be read.
-        #     particle_species_labels: N-BODY STAR
-        #     variable_labels: POSITION_X, VELOCITY_X, TIMESTEP
-        self.num_particle_species = param_handle.parameters['num_particle_species'][0]
-        self.num_primary_variables = param_handle.parameters['num_primary_variables']
-        self.num_secondary_variables = param_handle.parameters['num_secondary_variables']
-        # sfc values should come from subset and describe the domain for chunking
-        self.species_min = 0   
-        self.species_max = self.num_particle_species
- 
-        self.particle_variable_labels = {}
-        # combine primary and secondary varibles into one particle_variable labels 
-        # organized by species count
-        self.num_particle_variable_labels=[]
-        for ispec in xrange(self.num_particle_species) :
-            speclabelprimary = "species_%02d_primary_variable_labels" % ispec           
-            self.particle_variable_labels[ispec] = param_handle.parameters[speclabelprimary]
-            if self.num_secondary_variables[ispec] > 0 :
-                speclabelsecondary = "species_%02d_secondary_variable_labels" % ispec           
-                self.particle_variable_labels[ispec] = param_handle.parameters[speclabelprimary]+\
-                    param_handle.parameters[speclabelsecondary]
-                #self.num_particle_variable_labels[ispec]=self.num_primary_variables[ispec]+self.num_secondary_variables[ispec]
-                
-        # ####### variable dictionary from artio to yt ############
-        self.arttoyt_label_spec = {}
-        self.particle_species_labels = param_handle.parameters['particle_species_labels']
-        for label in self.particle_species_labels : 
-            print 'particle labels in caller.pyx',label
-            if label == 'N-BODY' : 
-                self.arttoyt_label_spec[label] = 'nbody' 
-            elif label == 'STARS' : 
-                self.arttoyt_label_spec[label] = 'stars' 
-            else :
-                self.arttoyt_label_spec[label] = label
-        self.arttoyt_label_var = {}
-        # 
-        for ispec in xrange(self.num_particle_species):
-            for label in self.particle_variable_labels[ispec] :
-                if label == 'BIRTH_TIME' :
-                    self.arttoyt_label_var[(ispec,label)] = 'creation_time'
-                else :
-                    self.arttoyt_label_var[(ispec,label)] = label
-                    print (ispec,label),label
-        print '----------------- arttoyt_label_spec (caller.pyx)'
-        print self.arttoyt_label_spec
-        print '----------------- arttoyt_label_var (caller.pyx)'
-        print self.arttoyt_label_var
-        print '-----------------'
-        print 'particle arttoyt translation'
-
-    def particle_pos_fill(self,accessed_species, fieldnames) :
-        #snl FIX the following should go in fields or be a subclass
-        #what are the matched art species? 
-        art_isubspec={}
-        art_labelispec={}
-        self.art_ispec={}
-        ispec=-1 
-        for artlabelspec in self.particle_species_labels : #N-BODY, N-BODY, STAR
-            ispec+=1
-            for ytlabelspec in accessed_species : 
-                if self.arttoyt_label_spec[artlabelspec] == ytlabelspec : 
-                    print 'match, in species labels?', ytlabelspec, accessed_species
-                    if art_isubspec[ytlabelspec] == None :
-                        isubspec = 0
-                    else : 
-                        isubspec += 1
-                    self.art_ispec.append(ispec)
-                    art_isubspec[ytlabelspec].append(isubspec)
-                    art_labelispec[ytlabelspec].append(ispec)
-#                    self.matched_labelspec.append(ytlabelspec)
-                    
-        # FIX Hardcoded subspecies selection
-        # e.g., only read least massive N-BODY particles -- pop off the species corresponding 
-        # to the subspecies you don't want
-        min_mass = 1e30
-        imin_mass = 0
-        for isubspec in art_isubspec['nbody'] : 
-            ispec = art_labelispec['nbody'][isubspec]
-            subspecmass = self.param_handle.parameters['particle_species_mass'][ispec] 
-            if subspecmass < min_mass :
-                min_mass = subspecmass
-                imin_mass = isubspec
-        # pop elements in nbody that are not min_mass
-        for isubspec in art_isubspec['nbody'] : 
-            ispec = art_labelispec['nbody'][isubspec]
-            if isubspec != imin_mass :
-                art_isubspec['nbody'].pop(ispec)
-                self.art_ispec.pop(ispec)
-                art_labelispec['nbody'].pop(ispec)
-        print 'art_isubspec[nbody] after pop...', art_isubspec['nbody']
-        print self.art_ispec
-        print art_labelispec['nbody'].pop(ispec)
-
-        print 'exiting in particle_pos_fill' 
-        sys.exit(1)
-         
-        #what are the matched art variable indices and (yt) labels for the matched species?
-        self.art_ivar=[]
-        self.matched_labelvar=[]
-        for ispec in self.art_ispec :
-            self.matched_labelvar[ispec]={}
-            i=-1
-            artivar=-1
-            for artlabelvar in self.particle_variable_labels[ispec] :  
-            #POSITION_X,... BIRTH_TIME, INITIAL_MASS, MASS, METALLICITY_SNII, METALLICITY_SNIa
-                artivar+=1
-                i+=1
-                for ytlabelvar in fieldnames : #fieldnames could be a function of species
-                    if self.arttoyt_label_var[artlabelvar] == ytlabelvar :
-                        print 'match, in fieldnames?',ispec, ytlabelvar, fieldnames
-                        self.art_ivar[(ispec,ytlabelvar)]=i
-#                        self.art_ivar[(ispec,artivar)]=i
-                        self.matched_labelvar[ispec].append(ytlabelvar)
-            print 'ispec=',ispec,'matched variable labels', self.matched_labelvar[ispec]
-            print 'art index of matched fieldnames',self.art_ivar
-        print 'exiting in particle_pos_fill' 
-        sys.exit(1)
-################ The above should be part of a subclass that uses accessed_species #######       
-
-        cdef artio_file handle
-        self.accessed_species = accessed_species
-        
-        self.selection={}
-        #ouch we are reading the particle file nspecies times?!
-        for ispec in self.art_ispec : 
-#              ispec->      self.species_min, self.species_max,\ 
-            handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                         ARTIO_OPEN_PARTICLES, artio_context_global ) 
-            status = artio_particle_read_sfc_range_yt(handle,\
-                    self.sfc_min, self.sfc_max,\
-                    ispec, ispec,\
-                    wrap_particle_pos_callback, <void*>self) 
-            check_artio_status(status, artio_grid_routines.__name__)
-            artio_fileset_close(handle) 
-            
-        print 'done reading particle positions'
-        return self.selection
-    def particle_var_fill(self, source, fieldnames, accessed_species, particle_mask):
-        cdef artio_file handle
-        self.source = source
-        if len(self.art_ivar) > 0 :
-            self.count=0
-            for ispec in self.art_ispec : 
-                handle = artio_fileset_open( self.param_handle.file_prefix, 
-                                             ARTIO_OPEN_PARTICLES, artio_context_global ) 
-                #snl check this:
-                status = artio_particle_read_sfc_range_yt(handle,\
-                        self.sfc_min, self.sfc_max,\
-                        self.species_min, self.species_max,\
-                        wrap_particle_var_callback, <void*>self)
-                check_artio_status(status, artio_grid_routines.__name__)
-                artio_fileset_close(handle) 
-        print 'done buffering variables'
-    def particle_pos_callback(self, pos):
-        self.selection['x'].append(pos[0])
-        self.selection['y'].append(pos[1])
-        self.selection['z'].append(pos[2])
-    def particle_var_callback(self, particle_var, ispec):
-        for labelvar in self.matched_labelvar[ispec] : 
-            artivar = self.art_ivar[(ispec,labelvar)]
-            self.source[labelvar][self.count] = particle_var[artivar] 
-        self.count=self.count+1
- 
-#        callback(pid,
-#                 primary_variables,
-#                 secondary_variables,
-#                 species, subspecies, sfc)
-###### callbacks (e.g. https://github.com/cython/cython/tree/master/Demos/callback) ######
-cdef void wrap_particle_pos_callback(int pid, double *primary_variables, 
-                                     double *secondary_variables, int species, 
-                                     int subspecies, int isfc, void *pyobject ) :
-    artioroutines = <object>pyobject
-    #FIX hardcoded assumption that position always comes first in primary variables
-    #assert this elsewhere
-    pos=[]
-    pos[0:2]=[primary_variables[0],primary_variables[1],primary_variables[2]]
-    artioroutines.particle_pos_callback(pos)
-
-cdef void wrap_particle_var_callback(int pid, double *primary_variables, double *secondary_variables, int species, 
-                                 int subspecies, int isfc, void *pyobject ) :
-    artioroutines = <object>pyobject
-    cell_var={}
-    # subspecies is only relevant when there are multiple "types" of "star-particles"  (e.g. AGN) 
-    ispec = species
-    particle_var = [ primary_variables[i] for i in range(artioroutines.num_primary_variables[ispec]) ]\
-        +        [ secondary_variables[i] for i in range(artioroutines.num_secondary_variables[ispec]) ]
-    print 'snl count particle var',(artioroutines.num_primary_variables[ispec]+artioroutines.num_secondary_variables[ispec]), particle_var 
-    artioroutines.cell_var_callback(particle_var)
-        
 ###################################################
 def artio_is_valid( char *file_prefix ) :
-    cdef artio_file handle = artio_fileset_open( file_prefix, 
+    cdef artio_fileset_handle *handle = artio_fileset_open( file_prefix, 
             ARTIO_OPEN_HEADER, artio_context_global )
     if handle == NULL :
         return False
     else :
         artio_fileset_close(handle) 
     return True
-

File yt/frontends/artio/artio_headers/artio.c

 #include "artio.h"
 #include "artio_internal.h"
 
-artio_file artio_fileset_allocate( char *file_prefix, int mode,
-        artio_context context );
-void artio_fileset_destroy( artio_file handle );
+artio_fileset *artio_fileset_allocate( char *file_prefix, int mode,
+		const artio_context *context );
+void artio_fileset_destroy( artio_fileset *handle );
 
-artio_file artio_fileset_open(char * file_prefix, int type, artio_context context) {
-	artio_fh head_fh;
+artio_fileset *artio_fileset_open(char * file_prefix, int type, const artio_context *context) {
+	artio_fh *head_fh;
 	char filename[256];
 	int ret;
+	int64_t tmp;
 
-	artio_file handle = 
+	artio_fileset *handle = 
 		artio_fileset_allocate( file_prefix, ARTIO_FILESET_READ, context );
 	if ( handle == NULL ) {
 		return NULL;
 		return NULL;
 	}
 
-	ret = artio_parameter_read(head_fh, &handle->param_list);
+	ret = artio_parameter_read(head_fh, handle->parameters );
 	if ( ret != ARTIO_SUCCESS ) {
 		artio_fileset_destroy(handle);
 		return NULL;
 	artio_file_fclose(head_fh);
 
 	artio_parameter_get_long(handle, "num_root_cells", &handle->num_root_cells);
+	
+	if ( artio_parameter_get_int(handle, "sfc_type", &handle->sfc_type ) != ARTIO_SUCCESS ) {
+		handle->sfc_type = ARTIO_SFC_HILBERT;
+	}
+
+	handle->nBitsPerDim = 0;
+	tmp = handle->num_root_cells >> 3;
+	while ( tmp ) {
+		handle->nBitsPerDim++;
+		tmp >>= 3;
+	}
 
 	/* default to accessing all sfc indices */
 	handle->proc_sfc_begin = 0;
 	return handle;
 }
 
-artio_file artio_fileset_create(char * file_prefix, int64_t root_cells, 
-		int64_t proc_sfc_begin, int64_t proc_sfc_end, artio_context context) {
-    artio_file handle = 
+artio_fileset *artio_fileset_create(char * file_prefix, int64_t root_cells, 
+		int64_t proc_sfc_begin, int64_t proc_sfc_end, const artio_context *context) {
+    artio_fileset *handle = 
 		artio_fileset_allocate( file_prefix, ARTIO_FILESET_WRITE, context );
     if ( handle == NULL ) {
         return NULL;
 
 #ifdef ARTIO_MPI
 	MPI_Allgather( &proc_sfc_begin, 1, MPI_LONG_LONG, 
-			handle->proc_sfc_index, 1, MPI_LONG_LONG, context->comm );
+			handle->proc_sfc_index, 1, MPI_LONG_LONG, handle->context->comm );
 #else
 	handle->proc_sfc_index[0] = 0;
 #endif /* ARTIO_MPI */
 	return handle;
 }
 
-int artio_fileset_close(artio_file handle) {
+int artio_fileset_close(artio_fileset *handle) {
 	char header_filename[256];
-	artio_fh head_fh;
+	artio_fh *head_fh;
 	
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 		}
 
 		if (0 == handle->rank) {
-			artio_parameter_write(head_fh, &handle->param_list);
+			artio_parameter_write(head_fh, handle->parameters );
 		}
 
 		artio_file_fclose(head_fh);
 	return ARTIO_SUCCESS;
 }
 
-artio_file artio_fileset_allocate( char *file_prefix, int mode, 
-		artio_context context ) {
+artio_fileset *artio_fileset_allocate( char *file_prefix, int mode, 
+		const artio_context *context ) {
 	int my_rank;
 	int num_procs;
 
-    artio_file handle = (artio_file)malloc(sizeof(struct artio_file_struct));
+    artio_fileset *handle = (artio_fileset *)malloc(sizeof(artio_fileset));
 	if ( handle != NULL ) {
-		artio_parameter_list_init(&handle->param_list);
+		handle->parameters = artio_parameter_list_init();
+
+		handle->context = (artio_context *)malloc(sizeof(artio_context));
+		if ( handle->context == NULL ) {
+			return NULL;
+		}
+		memcpy( handle->context, context, sizeof(artio_context) );
 
 #ifdef ARTIO_MPI
-		MPI_Comm_size(context->comm, &num_procs);
-		MPI_Comm_rank(context->comm, &my_rank);
+		MPI_Comm_size(handle->context->comm, &num_procs);
+		MPI_Comm_rank(handle->context->comm, &my_rank);
 #else
 		num_procs = 1;
 		my_rank = 0;
 
 		handle->rank = my_rank;
 		handle->num_procs = num_procs;
-		handle->context = context;
 		handle->endian_swap = 0;
 
 		handle->proc_sfc_index = NULL;
 	return handle;
 }
 
-void artio_fileset_destroy( artio_file handle ) {
+void artio_fileset_destroy( artio_fileset *handle ) {
 	if ( handle == NULL ) return;
 
 	if ( handle->proc_sfc_index != NULL ) free( handle->proc_sfc_index );
         artio_fileset_close_particles(handle);
     }
 
-	artio_parameter_free_list(&handle->param_list);
+	if ( handle->context != NULL ) free( handle->context );
+
+	artio_parameter_list_free(handle->parameters);
 
 	free(handle);
 }

File yt/frontends/artio/artio_headers/artio.h

 #define ARTIO_READ_LEAFS                    1
 #define ARTIO_READ_REFINED                  2
 #define	ARTIO_READ_ALL                      3
-#define	ARTIO_READ_REFINED_NOT_ROOT         4
+
+#define ARTIO_RETURN_OCTS					4
+#define ARTIO_RETURN_CELLS					0
 
 /* allocation strategy */
 #define ARTIO_ALLOC_EQUAL_SFC               0
 #define ARTIO_ALLOC_EQUAL_PROC              1
-#define ARTIO_ALLOC_MAX_FILE_SIZE  	    2
+#define ARTIO_ALLOC_MAX_FILE_SIZE  	        2
 
 #define ARTIO_TYPE_STRING                   0
 #define ARTIO_TYPE_CHAR                     1
 #define ARTIO_ERR_INVALID_SPECIES           108
 #define ARTIO_ERR_INVALID_ALLOC_STRATEGY    109
 #define ARTIO_ERR_INVALID_LEVEL             110
-#define ARTIO_ERR_INVALID_PARAM_LIST        111
+#define ARTIO_ERR_INVALID_PARAMETER_LIST    111
 #define ARTIO_ERR_INVALID_DATATYPE          112
 #define ARTIO_ERR_INVALID_OCT_REFINED       113
 #define ARTIO_ERR_INVALID_HANDLE            114
+#define ARTIO_ERR_INVALID_CELL_TYPES        115
 
 #define ARTIO_ERR_DATA_EXISTS               200
 #define ARTIO_ERR_INSUFFICIENT_DATA         201
 #ifdef ARTIO_MPI
 #include <mpi.h>
 
-struct artio_context_struct {
+typedef struct {
     MPI_Comm comm;
-};
+} artio_context;
 #else
-struct artio_context_struct {
+typedef struct {
     int comm;
-};
+} artio_context;
 #endif
 
-typedef struct artio_file_struct * artio_file;
-typedef struct artio_param_list * artio_parameters;
-typedef struct artio_context_struct * artio_context;
+typedef struct artio_fileset_struct artio_fileset;
 
-extern artio_context artio_context_global;
+extern const artio_context *artio_context_global;
 
 /*
  * Description: Open the file
  *
- *  filename			The file prefix
+ *  filename		The file prefix
  *  type			combination of ARTIO_OPEN_PARTICLES and ARTIO_OPEN_GRID flags
  */
-artio_file artio_fileset_open(char * file_name, int type, artio_context context);
+artio_fileset *artio_fileset_open( char * file_name, int type, const artio_context *context);
 
 /**
  * Description: Create fileset and begin populating header information
  *  handle			the artio file handle
  *
  */
-artio_file artio_fileset_create(char * file_prefix, 
-        int64_t root_cells, int64_t proc_sfc_begin, int64_t proc_sfc_end, artio_context context);
+artio_fileset *artio_fileset_create(char * file_prefix, 
+        int64_t root_cells, int64_t proc_sfc_begin, int64_t proc_sfc_end, const artio_context *context);
 
 /*
  * Description	Close the file
  */
-int artio_fileset_close(artio_file handle);
+int artio_fileset_close(artio_fileset *handle);
 
 /* public parameter interface */
-int artio_parameter_iterate( artio_file handle, char *key, int *type, int *length );
-int artio_parameter_get_array_length(artio_file handle, char * key, int *length);
+int artio_parameter_iterate( artio_fileset *handle, char *key, int *type, int *length );
+int artio_parameter_get_array_length(artio_fileset *handle, const char * key, int *length);
 
-void artio_parameter_set_int(artio_file handle, char * key, int32_t value);
-int artio_parameter_get_int(artio_file handle, char * key, int32_t * value);
+int artio_parameter_set_int(artio_fileset *handle, const char * key, int32_t value);
+int artio_parameter_get_int(artio_fileset *handle, const char * key, int32_t * value);
 
-void artio_parameter_set_int_array(artio_file handle, char * key, int length,
+int artio_parameter_set_int_array(artio_fileset *handle, const char * key, int length,
 		int32_t *values);
-int artio_parameter_get_int_array(artio_file handle, char * key, int length,
+int artio_parameter_get_int_array(artio_fileset *handle, const char * key, int length,
 		int32_t *values);
 
-void artio_parameter_set_string(artio_file handle, char * key, char * value);
-int artio_parameter_get_string(artio_file handle, char * key, char * value, int max_length);
+int artio_parameter_set_string(artio_fileset *handle, const char * key, char * value);
+int artio_parameter_get_string(artio_fileset *handle, const char * key, char * value, int max_length);
 
-void artio_parameter_set_string_array(artio_file handle, char * key,
+int artio_parameter_set_string_array(artio_fileset *handle, const char * key,
 		int length, char ** values);
-int artio_parameter_get_string_array(artio_file handle, char * key,
+int artio_parameter_get_string_array(artio_fileset *handle, const char * key,
 		int length, char ** values, int max_length);
 
-void artio_parameter_set_float(artio_file handle, char * key, float value);
-int artio_parameter_get_float(artio_file handle, char * key, float * value);
+int artio_parameter_set_float(artio_fileset *handle, const char * key, float value);
+int artio_parameter_get_float(artio_fileset *handle, const char * key, float * value);
 
-void artio_parameter_set_float_array(artio_file handle, char * key,
+int artio_parameter_set_float_array(artio_fileset *handle, const char * key,
 		int length, float *values);
-int artio_parameter_get_float_array(artio_file handle, char * key,
+int artio_parameter_get_float_array(artio_fileset *handle, const char * key,
 		int length, float * values);
 
-void artio_parameter_set_double(artio_file handle, char * key, double value);
-int  artio_parameter_get_double(artio_file handle, char * key, double * value);
+int artio_parameter_set_double(artio_fileset *handle, const char * key, double value);
+int  artio_parameter_get_double(artio_fileset *handle, const char * key, double * value);
 
-void artio_parameter_set_double_array(artio_file handle, char * key,
+int artio_parameter_set_double_array(artio_fileset *handle, const char * key,
 		int length, double * values);
-int artio_parameter_get_double_array(artio_file handle, char * key,
+int artio_parameter_get_double_array(artio_fileset *handle, const char * key,
         int length, double *values);
 
-void artio_parameter_set_long(artio_file handle, char * key, int64_t value);
-int artio_parameter_get_long(artio_file handle, char * key, int64_t *value);
+int artio_parameter_set_long(artio_fileset *handle, const char * key, int64_t value);
+int artio_parameter_get_long(artio_fileset *handle, const char * key, int64_t *value);
 
-void artio_parameter_set_long_array(artio_file handle, char * key,
+int artio_parameter_set_long_array(artio_fileset *handle, const char * key,
         int length, int64_t *values);
-int artio_parameter_get_long_array(artio_file handle, char * key,
+int artio_parameter_get_long_array(artio_fileset *handle, const char * key,
         int length, int64_t *values);
 
 /* public grid interface */
-typedef void (* GridCallBack)(float * variables, int level, int refined,
-		int64_t sfc_index);
-typedef void (* GridCallBackYTPos)(double * variables, int level, int refined,
-                                 int64_t sfc_index, double pos[3], void * pyobject);
-typedef void (* GridCallBackYT)(double * variables, int level, int refined,
-                                 int64_t sfc_index, void * pyobject);
+typedef void (* GridCallBack)( int64_t sfc_index, int level,
+		double *pos, float * variables, int *refined );
 
 /*
  * Description:	Add a grid component to a fileset open for writing
  *  num_levels_per_root_tree	Maximum tree depth for each oct tree
  *  num_octs_per_root_tree	Total octs in each oct tree
  */
-int artio_fileset_add_grid(artio_file handle,
+int artio_fileset_add_grid(artio_fileset *handle,
         int num_grid_files, int allocation_strategy,
         int num_grid_variables,
         char ** grid_variable_labels,
         int * num_levels_per_root_tree,
         int * num_octs_per_root_tree );
 
-int artio_fileset_open_grid(artio_file handle);
-int artio_fileset_close_grid(artio_file handle);
-int artio_fileset_open_particle(artio_file handle);
-int artio_fileset_close_particle(artio_file handle);
+int artio_fileset_open_grid(artio_fileset *handle);
+int artio_fileset_close_grid(artio_fileset *handle);
+int artio_fileset_open_particle(artio_fileset *handle);
+int artio_fileset_close_particle(artio_fileset *handle);
 
 /*
  * Description:	Output the variables of the root level cell and the hierarchy of the Oct tree correlated with this root level cell
  *  level			The depth of the Oct tree correlated to the root level cell
  *  num_level_octs		The array store the number of Oct nodes each level
  */
-int artio_grid_write_root_cell_begin(artio_file handle, int64_t sfc, 
+int artio_grid_write_root_cell_begin(artio_fileset *handle, int64_t sfc, 
 		float * variables, int level, int * num_octs_per_level);
 
 /*
  * Description:	Do something at the end of writing the root level cell
  */
-int artio_grid_write_root_cell_end(artio_file handle);
+int artio_grid_write_root_cell_end(artio_fileset *handle);
 
 /*
  * Description:	Do something at the beginning of each level
  */
-int artio_grid_write_level_begin(artio_file handle, int level );
+int artio_grid_write_level_begin(artio_fileset *handle, int level );
 
 /*
  * Description:	Do something at the end of each level
  */
-int artio_grid_write_level_end(artio_file handle);
+int artio_grid_write_level_end(artio_fileset *handle);
 
 /*
  * Description:	Output the data of a special oct tree node to the file
  *  handle			The handle of the file
  *  variables 			The array recording the variables of the eight cells belonging to this Octree node.
  */
-int artio_grid_write_oct(artio_file handle, float *variables, int *refined);
+int artio_grid_write_oct(artio_fileset *handle, float *variables, int *refined);
 
 /*
  * Description:	Read the variables of the root level cell and the hierarchy of the Octtree
  *  num_octs_per_level		The number of node of each oct level
  *
  */
-int artio_grid_read_root_nocts(artio_file handle, int64_t sfc, float *variables,
-		int *num_tree_levels, int *num_octs_per_level);
-int artio_grid_read_root_cell_begin(artio_file handle, int64_t sfc, float *variables,
+int artio_grid_read_root_cell_begin(artio_fileset *handle, int64_t sfc, 
+		double *pos, float *variables,
 		int *num_tree_levels, int *num_octs_per_level);
 
 /*
  * Description:	Do something at the end of reading the root level cell
  */
-int artio_grid_read_root_cell_end(artio_file handle);
+int artio_grid_read_root_cell_end(artio_fileset *handle);
 
 /*
  * Description:	Do something at the beginning of each level
  */
-int artio_grid_read_level_begin(artio_file handle, int level );
+int artio_grid_read_level_begin(artio_fileset *handle, int level );
 
 /*
  * Description:	Do something at the end of each level
  */
-int artio_grid_read_level_end(artio_file handle);
+int artio_grid_read_level_end(artio_fileset *handle);
 
 /*
  * Description:	Read the data of a special oct tree node from the file
  */
-int artio_grid_read_oct(artio_file handle, float *variables, int *refined);
+int artio_grid_read_oct(artio_fileset *handle, double *pos, 
+		float *variables, int *refined);
 
-int artio_grid_cache_sfc_range(artio_file handle, int64_t sfc_start, int64_t sfc_end);
+int artio_grid_cache_sfc_range(artio_fileset *handle, int64_t sfc_start, int64_t sfc_end);
+int artio_grid_clear_sfc_cache(artio_fileset *handle );
+
+int artio_grid_count_octs_in_sfc_range(artio_fileset *handle,
+        int64_t start, int64_t end, int64_t *num_octs);
 
 /*
  * Description:	Read a segment of oct nodes
  *  option			1. refined nodes; 2 leaf nodes; 3 all nodes
  *  callback			callback function
  */
-int artio_grid_read_sfc_range(       artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read, int max_level_to_read, int options, GridCallBack callback);
-int artio_grid_read_sfc_range_ytpos(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read,int max_level_to_read, int options, GridCallBackYTPos callback, void *pyobject);
-int artio_grid_read_sfc_range_yt(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read, int max_level_to_read, int options, GridCallBackYT callback, void *pyobject);
-		
+int artio_grid_read_sfc_range(artio_fileset *handle, int64_t sfc1, int64_t sfc2, 
+		int min_level_to_read, int max_level_to_read, 
+		int options, GridCallBack callback);
 
-typedef void (* ParticleCallBack)(int64_t pid, 
-		double *primary_variables, float *secondary_variables, 
-		int species, int subspecies, int64_t sfc_index);
-typedef void (* ParticleCallBackYT)(int64_t pid, 
-		double *primary_variables, float *secondary_variables, 
-                int species, int subspecies, int64_t sfc_index, void *pyobject);
+
+typedef void (* ParticleCallBack)(int64_t sfc_index,
+		int species, int subspecies, int64_t pid, 
+		double *primary_variables, float *secondary_variables );
 
 /**
  *  header			head file name
  *  handle			the artio file handle
  *
  */
-int artio_fileset_add_particles(artio_file handle, 
+int artio_fileset_add_particles(artio_fileset *handle, 
         int num_particle_files, int allocation_strategy,
         int num_species, char **species_labels,
         int *num_primary_variables,
         char ***secondary_variable_labels_per_species,
         int *num_particles_per_species_per_root_tree );
 
-int artio_fileset_open_particles(artio_file handle);
-int artio_fileset_close_particles(artio_file handle);
+int artio_fileset_open_particles(artio_fileset *handle);
+int artio_fileset_close_particles(artio_fileset *handle);
 
 /*
  * Description:	Output the variables of the root level cell and the hierarchy of 
  *  level			The depth of the Oct tree correlated to the root level cell
  *  num_level_octs		The array store the number of Oct nodes each level
  */
-int artio_particle_write_root_cell_begin(artio_file handle, int64_t sfc,
+int artio_particle_write_root_cell_begin(artio_fileset *handle, int64_t sfc,
 		int *num_particles_per_species);
 
 /*
  * Description:	Do something at the end of writing the root level cell
  */
-int artio_particle_write_root_cell_end(artio_file handle);
+int artio_particle_write_root_cell_end(artio_fileset *handle);
 
 /*
  * Description:	Do something at the beginning of each level
  */
-int artio_particle_write_species_begin(artio_file handle, int species );
+int artio_particle_write_species_begin(artio_fileset *handle, int species );
 
 /*
  * Description:	Do something at the end of each level
  */
-int artio_particle_write_species_end(artio_file handle);
+int artio_particle_write_species_end(artio_fileset *handle);
 
 /*
  * Description: Output the data of a special oct tree node to the file
  *  handle			The handle of the file
  *  variables 			The array recording the variables of the eight cells belonging to this Octree node.
  */
-int artio_particle_write_particle(artio_file handle, int64_t pid, int subspecies, 
+int artio_particle_write_particle(artio_fileset *handle, int64_t pid, int subspecies, 
 			double* primary_variables, float *secondary_variables);
 
 /*
  *  num_octs_per_level		The number of node of each oct level
  *
  */
-int artio_particle_read_root_cell_begin(artio_file handle, int64_t sfc, 
+int artio_particle_read_root_cell_begin(artio_fileset *handle, int64_t sfc, 
 			int * num_particle_per_species);
 
 /*
  * Description:	Do something at the end of reading the root level cell
  */
-int artio_particle_read_root_cell_end(artio_file handle);
+int artio_particle_read_root_cell_end(artio_fileset *handle);
 
 /*
  * Description:	Do something at the beginning of each level
  */
-int artio_particle_read_species_begin(artio_file handle, int species );
+int artio_particle_read_species_begin(artio_fileset *handle, int species );
 
 /*
  * Description:  Do something at the end of each level
  */
-int artio_particle_read_species_end(artio_file handle);
+int artio_particle_read_species_end(artio_fileset *handle);
 
 /*
  * Description:	Read the data of a single particle from the file
  */
-int artio_particle_read_particle(artio_file handle, int64_t *pid, int *subspecies,
+int artio_particle_read_particle(artio_fileset *handle, int64_t *pid, int *subspecies,
 			double *primary_variables, float *secondary_variables);
 
-int artio_particle_cache_sfc_range(artio_file handle, int64_t sfc_start, int64_t sfc_end);
+int artio_particle_cache_sfc_range(artio_fileset *handle, int64_t sfc_start, int64_t sfc_end);
 
 /*
  * Description: Read a segment of particles
  *  end_species			the last particle species to read
  *  callback			callback function
  */
-int artio_particle_read_sfc_range(artio_file handle, 
+int artio_particle_read_sfc_range(artio_fileset *handle, 
 		int64_t sfc1, int64_t sfc2, 
 		int start_species, int end_species,
 		ParticleCallBack callback);
 
-int artio_particle_read_sfc_range_yt(artio_file handle, 
-		int64_t sfc1, int64_t sfc2, 
-		int start_species, int end_species,
-                ParticleCallBackYT callback, void *pyobject);
-
 #endif /* __ARTIO_H__ */

File yt/frontends/artio/artio_headers/artio_grid.c

 #include <stdlib.h>
 #include <stdint.h>
 
-#include "sfc.h"
 #include "artio.h"
 #include "artio_internal.h"
 
-int artio_grid_find_file(artio_grid_file ghandle, int start, int end, int64_t sfc);
-artio_grid_file artio_grid_file_allocate(void);
-void artio_grid_file_destroy(artio_grid_file ghandle);
+int artio_grid_find_file(artio_grid_file *ghandle, int start, int end, int64_t sfc);
+artio_grid_file *artio_grid_file_allocate(void);
+void artio_grid_file_destroy(artio_grid_file *ghandle);
+
+const double oct_pos_offsets[8][3] = {
+	{ -0.5, -0.5, -0.5 }, {  0.5, -0.5, -0.5 }, 
+	{ -0.5,  0.5, -0.5 }, {  0.5,  0.5, -0.5 }, 
+	{ -0.5, -0.5,  0.5 }, {  0.5, -0.5,  0.5 },
+	{ -0.5,  0.5,  0.5 }, {  0.5,  0.5,  0.5 }
+};
 
 /*
  * Open grid component of the fileset
  */
-int artio_fileset_open_grid(artio_file handle) {
+int artio_fileset_open_grid(artio_fileset *handle) {
 	int i;
 	char filename[256];
 	int first_file, last_file;
 	int mode;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 		return ARTIO_ERR_MEMORY_ALLOCATION;
 	}
 
-	ghandle->ffh = (artio_fh *)malloc(ghandle->num_grid_files * sizeof(artio_fh));
+	ghandle->ffh = (artio_fh **)malloc(ghandle->num_grid_files * sizeof(artio_fh *));
 	if ( ghandle->ffh == NULL ) {
 		artio_grid_file_destroy(ghandle);
 		return ARTIO_ERR_MEMORY_ALLOCATION;
 	return ARTIO_SUCCESS;
 }
 
-int artio_fileset_add_grid(artio_file handle, 
+int artio_fileset_add_grid(artio_fileset *handle, 
 		int num_grid_files, int allocation_strategy, 
 		int num_grid_variables, 
 		char ** grid_variable_labels,
 	char filename[256];
 	int mode;
 	int ret;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	}
 
 	/* allocate file handles */
-	ghandle->ffh = (artio_fh *)malloc(num_grid_files * sizeof(artio_fh));
+	ghandle->ffh = (artio_fh **)malloc(num_grid_files * sizeof(artio_fh *));
 	if ( ghandle->ffh == NULL ) {
 		artio_grid_file_destroy(ghandle);
 		return ARTIO_ERR_MEMORY_ALLOCATION;
 		}
 	}
 
-	handle->grid = ghandle;
+		handle->grid = ghandle;
 
 	artio_parameter_set_long_array(handle, "grid_file_sfc_index",
 			ghandle->num_grid_files + 1, ghandle->file_sfc_index);
 	return ARTIO_SUCCESS;
 }
 
-artio_grid_file artio_grid_file_allocate(void) {
-	artio_grid_file ghandle = 
-			(artio_grid_file)malloc(sizeof(struct artio_grid_file_struct));                                          
+artio_grid_file *artio_grid_file_allocate(void) {
+	artio_grid_file *ghandle = 
+			(artio_grid_file *)malloc(sizeof(struct artio_grid_file_struct));                                          
 	if ( ghandle != NULL ) {
 		ghandle->ffh = NULL;
 		ghandle->num_grid_variables = -1;
 		ghandle->cur_octs = -1;
 		ghandle->cur_sfc = -1;
 		ghandle->octs_per_level = NULL;
+
+		ghandle->pos_flag = 0;
+		ghandle->pos_cur_level = -1;
+		ghandle->next_level_size = -1;
+		ghandle->cur_level_size = -1;
+		ghandle->cell_size_level = 1e20;
+		ghandle->next_level_pos = NULL;
+		ghandle->cur_level_pos = NULL;
+		ghandle->next_level_oct = -1;
     }
 	return ghandle;
 }
 
-void artio_grid_file_destroy(artio_grid_file ghandle) {
+void artio_grid_file_destroy(artio_grid_file *ghandle) {
 	int i;
 	if ( ghandle == NULL ) return;	
 
 	if ( ghandle->sfc_offset_table != NULL ) free(ghandle->sfc_offset_table);
 	if ( ghandle->octs_per_level != NULL ) free(ghandle->octs_per_level);
 	if ( ghandle->file_sfc_index != NULL ) free(ghandle->file_sfc_index);
+	if ( ghandle->next_level_pos != NULL ) free(ghandle->next_level_pos);
+	if ( ghandle->cur_level_pos != NULL ) free(ghandle->cur_level_pos);
+
 	free(ghandle);
 }
 
-int artio_fileset_close_grid(artio_file handle) {
+int artio_fileset_close_grid(artio_fileset *handle) {
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	}
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_cache_sfc_range(artio_file handle, int64_t start, int64_t end) {
+int artio_grid_count_octs_in_sfc_range(artio_fileset *handle, 
+		int64_t start, int64_t end, int64_t *num_octs) {
+    int i;
+	int ret;
+	int file, first;
+	int64_t sfc;
+	int64_t offset, next_offset, size_offset;
+	int num_oct_levels;
+	int *num_octs_per_level;
+    artio_grid_file *ghandle;
+
+    if ( handle == NULL ) {
+        return ARTIO_ERR_INVALID_HANDLE;
+    }
+
+    if (handle->open_mode != ARTIO_FILESET_READ ||
+            !(handle->open_type & ARTIO_OPEN_GRID) ||
+            handle->grid == NULL ) {
+        return ARTIO_ERR_INVALID_FILESET_MODE;
+    }
+
+	if ( start > end || start < handle->proc_sfc_begin ||
+			end > handle->proc_sfc_end ) {
+		return ARTIO_ERR_INVALID_SFC_RANGE;
+	}
+
+	ghandle = handle->grid;
+
+	/* check that we're not in the middle of a read */
+	if ( ghandle->cur_sfc != -1 ) {
+		return ARTIO_ERR_INVALID_STATE;
+	}
+
+	*num_octs = 0;
+
+	if ( 8*ghandle->num_grid_variables <= ghandle->file_max_level ) {
+		/* we can't compute the number of octs through the offset table */
+		ret = artio_grid_cache_sfc_range( handle, start, end );
+		if ( ret != ARTIO_SUCCESS ) return ret;
+
+		num_octs_per_level = (int *)malloc(ghandle->file_max_level*sizeof(int) );
+		if ( num_octs_per_level == NULL ) {
+			return ARTIO_ERR_MEMORY_ALLOCATION;
+		}
+
+		for ( sfc = start; sfc <= end; sfc++ ) {
+			ret = artio_grid_read_root_cell_begin( handle, sfc, NULL, NULL, 
+					&num_oct_levels, num_octs_per_level );
+			if ( ret != ARTIO_SUCCESS ) return ret;
+
+			for ( i = 0; i < num_oct_levels; i++ ) {
+				*num_octs += num_octs_per_level[i];
+			}	
+
+			ret = artio_grid_read_root_cell_end( handle );
+			if ( ret != ARTIO_SUCCESS ) return ret;
+		}
+
+		free( num_octs_per_level );
+	} else {
+		/* TODO: add optimization if sfc range already cached */
+		file = artio_grid_find_file(ghandle, 0, ghandle->num_grid_files, start);
+		first = MAX( 0, start - ghandle->file_sfc_index[file] );
+
+		ret = artio_file_fseek(ghandle->ffh[file], 
+				sizeof(int64_t) * first, ARTIO_SEEK_SET);
+		if ( ret != ARTIO_SUCCESS ) return ret;
+
+		ret = artio_file_fread(ghandle->ffh[file], &offset, 1, ARTIO_TYPE_LONG );
+		if ( ret != ARTIO_SUCCESS ) return ret;
+
+		sfc = start;
+		while ( sfc <= end ) {
+		/* read next offset or compute end of file*/
+			if ( sfc < ghandle->file_sfc_index[file+1] - 1 ) {
+				ret = artio_file_fread(ghandle->ffh[file], &size_offset, 1, ARTIO_TYPE_LONG );
+				if ( ret != ARTIO_SUCCESS ) return ret;
+				next_offset = size_offset;
+			} else {
+				/* need to seek and ftell */
+				artio_file_fseek( ghandle->ffh[file], 0, ARTIO_SEEK_END );
+				artio_file_ftell( ghandle->ffh[file], &size_offset );
+				file++;
+
+				if ( sfc < end && file < ghandle->num_grid_files ) {
+					if ( sfc != ghandle->file_sfc_index[file] ) {
+						printf("ERROR!!!");
+						exit(1);
+					}
+					artio_file_fseek( ghandle->ffh[file], 0, ARTIO_SEEK_SET );
+					ret = artio_file_fread(ghandle->ffh[file], &next_offset, 
+							1, ARTIO_TYPE_LONG );
+					if ( ret != ARTIO_SUCCESS ) return ret;
+				}
+			}
+
+			/* this assumes (num_levels_per_root_tree)*sizeof(int) <
+			 *   size of an oct, or 8*num_variables > max_level so the 
+			 *   number of levels drops off in rounding to int */
+			*num_octs += (size_offset - offset - 
+				sizeof(float)*ghandle->num_grid_variables - sizeof(int) ) /
+				(8*(sizeof(float)*ghandle->num_grid_variables + sizeof(int) ));
+			offset = next_offset;
+			sfc++;
+		}
+	}
+
+	return ARTIO_SUCCESS;
+}
+
+int artio_grid_cache_sfc_range(artio_fileset *handle, int64_t start, int64_t end) {
 	int i;
 	int ret;
 	int first_file, last_file;
 	int64_t first, count, cur;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 		return ARTIO_ERR_INVALID_FILESET_MODE;
 	}
 
-	ghandle = handle->grid;
-
 	if ( start > end || start < handle->proc_sfc_begin ||
 			end > handle->proc_sfc_end ) {
 		return ARTIO_ERR_INVALID_SFC_RANGE;
 	}
 
+	artio_grid_clear_sfc_cache(handle);
 
-	if (ghandle->sfc_offset_table != NULL) {
-		free(ghandle->sfc_offset_table);
-	}
+	ghandle = handle->grid;
 
 	first_file = artio_grid_find_file(ghandle, 0, ghandle->num_grid_files, start);
 	last_file = artio_grid_find_file(ghandle, first_file, ghandle->num_grid_files, end);
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_find_file(artio_grid_file ghandle, int start, int end, int64_t sfc) {
+int artio_grid_clear_sfc_cache( artio_fileset *handle ) {
+    artio_grid_file *ghandle;
+
+    if ( handle == NULL ) {
+        return ARTIO_ERR_INVALID_HANDLE;
+    }
+
+    if (handle->open_mode != ARTIO_FILESET_READ ||
+            !(handle->open_type & ARTIO_OPEN_GRID) ||
+            handle->grid == NULL ) {
+        return ARTIO_ERR_INVALID_FILESET_MODE;
+    }
+
+    ghandle = handle->grid;
+
+	if ( ghandle->sfc_offset_table != NULL ) {
+		free(ghandle->sfc_offset_table);
+		ghandle->sfc_offset_table = NULL;
+	}
+
+    ghandle->cache_sfc_begin = -1;
+    ghandle->cache_sfc_end = -1;
+
+	return ARTIO_SUCCESS;
+}
+
+int artio_grid_find_file(artio_grid_file *ghandle, int start, int end, int64_t sfc) {
 	int j;
 
-	if ( ghandle == NULL ) {
-		return ARTIO_ERR_INVALID_HANDLE;
-	}
-
 	if ( start < 0 || start > ghandle->num_grid_files ||
 			end < 0 || end > ghandle->num_grid_files || 
 			sfc < ghandle->file_sfc_index[start] ||
 	}
 }
 
-int artio_grid_seek_to_sfc(artio_file handle, int64_t sfc) {
+int artio_grid_seek_to_sfc(artio_fileset *handle, int64_t sfc) {
 	int64_t offset;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 			offset, ARTIO_SEEK_SET);
 }
 
-int artio_grid_write_root_cell_begin(artio_file handle, int64_t sfc,
+int artio_grid_write_root_cell_begin(artio_fileset *handle, int64_t sfc,
 		float *variables, int num_oct_levels, int *num_octs_per_level) {
 	int i;
 	int ret;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_write_root_cell_end(artio_file handle) {
+int artio_grid_write_root_cell_end(artio_fileset *handle) {
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	}
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_write_level_begin(artio_file handle, int level) {
-	artio_grid_file ghandle;
+int artio_grid_write_level_begin(artio_fileset *handle, int level) {
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_write_level_end(artio_file handle) {
-	artio_grid_file ghandle;
+int artio_grid_write_level_end(artio_fileset *handle) {
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_write_oct(artio_file handle, float *variables,
+int artio_grid_write_oct(artio_fileset *handle, float *variables,
 		int *cellrefined) {
 	int i;
 	int ret;
-	artio_grid_file ghandle;
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	return ARTIO_SUCCESS;
 }
 
-int artio_grid_read_root_nocts(artio_file handle, int64_t sfc,
-		float *variables, int *num_oct_levels, int *num_octs_per_level) {
+/*
+ *
+ */
+int artio_grid_read_root_cell_begin(artio_fileset *handle, int64_t sfc,
+		double *pos, float *variables, int *num_oct_levels, 
+		int *num_octs_per_level) {
 	int i;
 	int ret;
-	artio_grid_file ghandle;
+	int coords[3];
+	artio_grid_file *ghandle;
 
 	if ( handle == NULL ) {
 		return ARTIO_ERR_INVALID_HANDLE;
 	ret = artio_grid_seek_to_sfc(handle, sfc);
 	if ( ret != ARTIO_SUCCESS ) return ret;
 
-	ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
-			variables, ghandle->num_grid_variables,
-			ARTIO_TYPE_FLOAT);
-	if ( ret != ARTIO_SUCCESS ) return ret;
+	if ( variables == NULL ) {
+		ret = artio_file_fseek( ghandle->ffh[ghandle->cur_file],
+				ghandle->num_grid_variables*sizeof(float), 
+				ARTIO_SEEK_CUR );
+		if ( ret != ARTIO_SUCCESS ) return ret;
+	} else {
+		ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
+				variables, ghandle->num_grid_variables,
+				ARTIO_TYPE_FLOAT);
+		if ( ret != ARTIO_SUCCESS ) return ret;
+	}
 
 	ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
 			num_oct_levels, 1, ARTIO_TYPE_INT);
 		return ARTIO_ERR_INVALID_OCT_LEVELS;
 	}
 
-	if (*num_oct_levels > 0) {
-		ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
-				num_octs_per_level, *num_oct_levels,
-				ARTIO_TYPE_INT);
-		if ( ret != ARTIO_SUCCESS ) return ret;
+	if ( pos != NULL ) {
+		ghandle->pos_flag = 1;
 
-		for (i = 0; i < *num_oct_levels; i++) {
-			ghandle->octs_per_level[i] = num_octs_per_level[i];
+		/* compute position from sfc */
+		artio_sfc_coords( handle, sfc, coords );
+		for ( i = 0; i < 3; i++ ) {
+			pos[i] = (double)coords[i] + 0.5;
 		}
-	}
 
-	ghandle->cur_sfc = sfc;
-	ghandle->cur_num_levels = *num_oct_levels;
-	ghandle->cur_level = -1;
+		if ( *num_oct_levels > 0 ) {
+			/* compute next level position */
+			if ( ghandle->next_level_pos == NULL ) {
+				ghandle->next_level_pos = (double *)malloc( 3*sizeof(double) );
+				if ( ghandle->next_level_pos == NULL ) {
+					return ARTIO_ERR_MEMORY_ALLOCATION;
+				}
+				ghandle->next_level_size = 1;
+			}
 
-	return ARTIO_SUCCESS;
-}
-int artio_grid_read_root_cell_begin(artio_file handle, int64_t sfc,
-		float *variables, int *num_oct_levels, int *num_octs_per_level) {
-	int i;
-	int ret;
-	artio_grid_file ghandle;
-
-	if ( handle == NULL ) {
-		return ARTIO_ERR_INVALID_HANDLE;
-	}
-
-	if (handle->open_mode != ARTIO_FILESET_READ || 
-			!(handle->open_type & ARTIO_OPEN_GRID) ||
-			handle->grid == NULL ) {
-		return ARTIO_ERR_INVALID_FILESET_MODE;
-	}
-
-	ghandle = handle->grid;
-
-	ret = artio_grid_seek_to_sfc(handle, sfc);
-	if ( ret != ARTIO_SUCCESS ) return ret;
-
-	ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
-			variables, ghandle->num_grid_variables,
-			ARTIO_TYPE_FLOAT);
-	if ( ret != ARTIO_SUCCESS ) return ret;
-
-	ret = artio_file_fread(ghandle->ffh[ghandle->cur_file], 
-			num_oct_levels, 1, ARTIO_TYPE_INT);
-	if ( ret != ARTIO_SUCCESS ) return ret;
-
-	if ( *num_oct_levels > ghandle->file_max_level ) {
-		return ARTIO_ERR_INVALID_OCT_LEVELS;
+			for ( i = 0; i < 3; i++ ) {
+				ghandle->next_level_pos[i] = pos[i];
+			}
+			ghandle->pos_cur_level = 0;
+		} else {
+			ghandle->pos_cur_level = -1;
+		}
+	} else {
+		ghandle->pos_flag = 0;
 	}
 
 	if (*num_oct_levels > 0) {
 	return ARTIO_SUCCESS;
 }
 
-/* Description  */
-int artio_grid_read_oct(artio_file handle, float *variables,
+int artio_grid_read_oct(artio_fileset *handle, 
+		double *pos,
+		float *variables,
 		int *refined) {
+	int i, j;
 	int ret;
-	artio_grid_file ghandle;
+	int local_refined[8];
+	artio_grid_file *ghandle;