Commits

Sam Leitner committed abb6c00

new artio frontend compiles with the rest of yt.... it is just a copy of ramses routines with some wrappers mixed in to _artio_reader.pyx for the artio libraries in artio_headers. A chunk of the oct routines were copied from yt/geometry as well.

Comments (0)

Files changed (25)

yt/frontends/artio/__config__.py

+# This file is generated by /Users/sleitner/repos/ytsnlupdate/setup.py
+# It contains system_info results at the time of building this package.
+__all__ = ["get_info","show"]
+
+
+def get_info(name):
+    g = globals()
+    return g.get(name, g.get(name + "_info", {}))
+
+def show():
+    for name,info_dict in globals().items():
+        if name[0] == "_" or type(info_dict) is not type({}): continue
+        print(name + ":")
+        if not info_dict:
+            print("  NOT AVAILABLE")
+        for k,v in info_dict.items():
+            v = str(v)
+            if k == "sources" and len(v) > 200:
+                v = v[:60] + " ...\n... " + v[-60:]
+            print("    %s = %s" % (k,v))
+    

yt/frontends/artio/__init__.py

 """
-API for yt.frontends.ramses
+API for yt.frontends.artio
 
 Author: Matthew Turk <matthewturk@gmail.com>
 Affiliation: UCSD

yt/frontends/artio/_artio_reader.pyx

   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-# Cython wrapping code for Oliver Hahn's RAMSES reader
+# Cython wrapping code for Oliver Hahn's ARTIO reader
 from cython.operator cimport dereference as deref, preincrement as inc
 from libc.stdlib cimport malloc, free, abs, calloc, labs
 
         char *c_str()
         string operator*()
 
-cdef extern from "RAMSES_typedefs.h":
+cdef extern from "ARTIO_typedefs.h":
     pass
 
-cdef extern from "RAMSES_info.hh" namespace "RAMSES":
+cdef extern from "ARTIO_info.hh" namespace "ARTIO":
     enum codeversion:
         version1
         version2
     #void hilbert3d(vector[double] &x, vector[double] &y, vector[double] &z,
     #               vector[double] &order, unsigned bit_length)
 
-cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::AMR":
+cdef extern from "ARTIO_amr_data.hh" namespace "ARTIO::AMR":
     cdef cppclass vec[real_t]:
         real_t x, y, z
         vec( real_t x_, real_t y_, real_t z_)
 
         bint is_refined(int ison)
 
-    cdef cppclass RAMSES_cell:
+    cdef cppclass ARTIO_cell:
         unsigned m_neighbour[6]
         unsigned m_father
         unsigned m_son[8]
 
         char m_pos
 
-        RAMSES_cell()
+        ARTIO_cell()
 
         bint is_refined(int ison)
 
         Cell_& operator[]( unsigned i)
         unsigned size()
 
-    cdef cppclass RAMSES_level:
+    cdef cppclass ARTIO_level:
         unsigned m_ilevel
-        vector[RAMSES_cell] m_level_cells
+        vector[ARTIO_cell] m_level_cells
         double m_xc[8]
         double m_yc[8]
         double m_zc[8]
 
         double m_dx
         unsigned m_nx
-        RAMSES_level (unsigned ilevel)
+        ARTIO_level (unsigned ilevel)
 
-        void register_cell( RAMSES_cell cell )
-        vector[RAMSES_cell].iterator begin()
-        vector[RAMSES_cell].iterator end()
+        void register_cell( ARTIO_cell cell )
+        vector[ARTIO_cell].iterator begin()
+        vector[ARTIO_cell].iterator end()
 
-        RAMSES_cell& operator[]( unsigned i)
+        ARTIO_cell& operator[]( unsigned i)
         unsigned size()
 
     # This class definition is out of date.  I have unrolled the template
 
         tree (snapshot &snap, int cpu, int maxlevel, int minlevel = 1)
 
-    cppclass tree_iterator "RAMSES::AMR::RAMSES_tree::iterator":
+    cppclass tree_iterator "ARTIO::AMR::ARTIO_tree::iterator":
         tree_iterator operator*()
-        RAMSES_cell operator*()
+        ARTIO_cell operator*()
         tree_iterator begin()
         tree_iterator end()
         tree_iterator to_parent()
         int get_absolute_position()
         int get_domain()
 
-    cdef cppclass RAMSES_tree:
+    cdef cppclass ARTIO_tree:
         # This is, strictly speaking, not a header.  But, I believe it is
         # going to work alright.
         cppclass header:
             double t
             double mass_sph
 
-        vector[RAMSES_level] m_AMR_levels
+        vector[ARTIO_level] m_AMR_levels
         vector[unsigned] mheadl, m_numbl, m_taill
 
         int m_cpu
         # typedef proto_iterator<const tree*> const_iterator;
         # typedef proto_iterator<tree *> iterator;
 
-        RAMSES_tree(snapshot &snap, int cpu, int maxlevel, int minlevel)
+        ARTIO_tree(snapshot &snap, int cpu, int maxlevel, int minlevel)
         void read()
 
         tree_iterator begin(int ilevel)
         vec[float] cell_pos_float "cell_pos<float>" (tree_iterator it, unsigned ind) 
         vec[float] grid_pos_float "grid_pos<float>" (tree_iterator it)
 
-cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::HYDRO":
+cdef extern from "ARTIO_amr_data.hh" namespace "ARTIO::HYDRO":
     enum hydro_var:
         density
         velocity_x
         pressure
         metallicit
 
-    char ramses_hydro_variables[][64]
+    char artio_hydro_variables[][64]
 
     # There are a number of classes here that we could wrap and utilize.
     # However, I will only wrap the methods I know we need.
         #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
         void read(string varname)
 
-    cdef cppclass RAMSES_hydro_data:
+    cdef cppclass ARTIO_hydro_data:
         cppclass header:
             unsigned ncpu
             unsigned nvar
         vector[string] m_varnames
         map[string, unsigned] m_var_name_map
 
-        RAMSES_hydro_data(RAMSES_tree &AMRtree)
+        ARTIO_hydro_data(ARTIO_tree &AMRtree)
         #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
         void read(string varname)
 
         vector[vector[double]] m_var_array
 
-cdef class RAMSES_tree_proxy:
+cdef class ARTIO_tree_proxy:
     cdef string *snapshot_name
     cdef snapshot *rsnap
 
-    cdef RAMSES_tree** trees
-    cdef RAMSES_hydro_data*** hydro_datas
+    cdef ARTIO_tree** trees
+    cdef ARTIO_hydro_data*** hydro_datas
 
     cdef int **loaded
 
     
     def __cinit__(self, char *fn):
         cdef int idomain, ifield, ii
-        cdef RAMSES_tree *local_tree
-        cdef RAMSES_hydro_data *local_hydro_data
+        cdef ARTIO_tree *local_tree
+        cdef ARTIO_hydro_data *local_hydro_data
         self.snapshot_name = new string(fn)
         self.rsnap = new snapshot(deref(self.snapshot_name), version3)
         # We now have to get our field names to fill our array
-        self.trees = <RAMSES_tree**>\
-            malloc(sizeof(RAMSES_tree*) * self.rsnap.m_header.ncpu)
+        self.trees = <ARTIO_tree**>\
+            malloc(sizeof(ARTIO_tree*) * self.rsnap.m_header.ncpu)
         for ii in range(self.ndomains): self.trees[ii] = NULL
-        self.hydro_datas = <RAMSES_hydro_data ***>\
-                       malloc(sizeof(RAMSES_hydro_data**) * self.rsnap.m_header.ncpu)
+        self.hydro_datas = <ARTIO_hydro_data ***>\
+                       malloc(sizeof(ARTIO_hydro_data**) * self.rsnap.m_header.ncpu)
         self.ndomains = self.rsnap.m_header.ncpu
         
         # Note we don't do ncpu + 1
         for idomain in range(self.rsnap.m_header.ncpu):
             # we don't delete local_tree
-            local_tree = new RAMSES_tree(deref(self.rsnap), idomain + 1,
+            local_tree = new ARTIO_tree(deref(self.rsnap), idomain + 1,
                                          self.rsnap.m_header.levelmax, 0)
             local_tree.read()
-            local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
-            self.hydro_datas[idomain] = <RAMSES_hydro_data **>\
-                malloc(sizeof(RAMSES_hydro_data*) * local_hydro_data.m_nvars)
+            local_hydro_data = new ARTIO_hydro_data(deref(local_tree))
+            self.hydro_datas[idomain] = <ARTIO_hydro_data **>\
+                malloc(sizeof(ARTIO_hydro_data*) * local_hydro_data.m_nvars)
             for ii in range(local_hydro_data.m_nvars):
                 self.hydro_datas[idomain][ii] = \
-                    new RAMSES_hydro_data(deref(local_tree))
+                    new ARTIO_hydro_data(deref(local_tree))
             self.trees[idomain] = local_tree
             # We do not delete the final snapshot, which we'll use later
             #if idomain + 1 < self.rsnap.m_header.ncpu:
         cdef int idomain, ifield
         # To ensure that 'delete' is used, not 'free',
         # we allocate temporary variables.
-        cdef RAMSES_tree *temp_tree
-        cdef RAMSES_hydro_data *temp_hdata
+        cdef ARTIO_tree *temp_tree
+        cdef ARTIO_hydro_data *temp_hdata
         for idomain in range(self.ndomains):
             for ifield in range(self.nfields):
                 temp_hdata = self.hydro_datas[idomain][ifield]
         # We need to do simulation domains here
 
         cdef unsigned idomain, ilevel
-        cdef RAMSES_tree *local_tree
-        cdef RAMSES_hydro_data *local_hydro_data
-        cdef RAMSES_level *local_level
+        cdef ARTIO_tree *local_tree
+        cdef ARTIO_hydro_data *local_hydro_data
+        cdef ARTIO_level *local_level
         cdef tree_iterator grid_it, grid_end
 
         # All the loop-local pointers must be declared up here
         # We delete and re-create
         cdef int varindex = self.field_ind[varname]
         cdef string *field_name = new string(varname)
-        cdef RAMSES_hydro_data *temp_hdata
+        cdef ARTIO_hydro_data *temp_hdata
         if self.loaded[domain_index][varindex] == 1:
             temp_hdata = self.hydro_datas[domain_index][varindex]
             del temp_hdata
             self.hydro_datas[domain_index][varindex] = \
-                new RAMSES_hydro_data(deref(self.trees[domain_index]))
+                new ARTIO_hydro_data(deref(self.trees[domain_index]))
             self.loaded[domain_index][varindex] = 0
         del field_name
 
         cdef unsigned idomain, ilevel
 
         # All the loop-local pointers must be declared up here
-        cdef RAMSES_tree *local_tree
-        cdef RAMSES_hydro_data *local_hydro_data
+        cdef ARTIO_tree *local_tree
+        cdef ARTIO_hydro_data *local_hydro_data
         cdef unsigned father
 
         cdef tree_iterator grid_it, grid_end, father_it
         cdef tree_iterator grid_it, grid_end
         cdef double* data = <double*> tr.data
 
-        cdef RAMSES_tree *local_tree = self.trees[domain - 1]
-        cdef RAMSES_hydro_data *local_hydro_data = self.hydro_datas[domain - 1][varindex]
+        cdef ARTIO_tree *local_tree = self.trees[domain - 1]
+        cdef ARTIO_hydro_data *local_hydro_data = self.hydro_datas[domain - 1][varindex]
 
         #inline ValueType_& cell_value( const typename TreeType_::iterator& it,
         #                               int ind )
                   np.ndarray[np.int32_t, ndim=3] filled,
                   int level, int ref_factor,
                   np.ndarray[np.int64_t, ndim=2] component_grid_info):
-        cdef RAMSES_tree *local_tree = NULL
-        cdef RAMSES_hydro_data *local_hydro_data = NULL
+        cdef ARTIO_tree *local_tree = NULL
+        cdef ARTIO_hydro_data *local_hydro_data = NULL
 
         cdef int gi, i, j, k, domain, offset
         cdef int ir, jr, kr
     else:
         rv_r = [R]
     return rv_r + rv_l
+
+
+##############################snl snl ARTIO reader
+from cython.operator cimport dereference as deref, preincrement as inc
+from libc.stdlib cimport malloc, free, abs, calloc, labs
+
+import numpy as np
+#cimport numpy as np
+cimport cython
+
+cdef extern from "math.h":
+    double ceil(double x)
+
+#struct artio_context_struct {
+#        MPI_Comm comm;
+#};
+cdef struct artio_context_struct: #in artio_mpi.h MPI_Comm
+    int comm
+
+#typedef struct artio_context_struct * artio_context;
+ctypedef artio_context_struct * artio_context
+
+cdef extern from "artio_headers/artio.c":
+
+     ctypedef int int64_t 
+
+     ctypedef struct param: 
+         int key_length
+         char key[64]
+         int val_length
+         int type 
+         char *value
+         
+     ctypedef struct list: 
+         param * head
+         param * tail
+         param * cursor
+         int iterate_flag 
+         int endian_swap 
+
+     #    artio_internal.h
+     ctypedef struct artio_particle_file:
+        #	artio_fh *ffh
+        int num_particle_files
+        int64_t *file_sfc_index
+        int64_t cache_sfc_begin
+        int64_t cache_sfc_end
+        int64_t *sfc_offset_table
+        
+        #/* maintained for consistency and user-error detection */
+        int num_species
+        int cur_file
+        int cur_species
+        int cur_particle
+        int64_t cur_sfc
+        int *num_primary_variables
+        int *num_secondary_variables
+        int *num_particles_per_species
+        
+     ctypedef struct artio_grid_file:
+        #        artio_fh *ffh
+        int num_grid_variables
+        int num_grid_files
+        int64_t *file_sfc_index
+        int64_t cache_sfc_begin
+        int64_t cache_sfc_end
+        int64_t *sfc_offset_table
+        
+        int file_max_level
+        #/* maintained for consistency and user-error detection */
+        int cur_file
+        int cur_num_levels
+        int cur_level
+        int cur_octs
+        int64_t cur_sfc
+        int *octs_per_level
+        
+        
+     ctypedef struct artio_file:
+        char file_prefix[256]
+        int endian_swap
+        int open_type
+        int open_mode
+        int rank
+        int num_procs
+        artio_context context #artio_mpi.h MPI_Comm
+        
+        int64_t *proc_sfc_index
+        int64_t proc_sfc_begin
+        int64_t proc_sfc_end
+        int64_t num_root_cells
+        
+        list param_list
+        artio_grid_file grid
+        artio_particle_file particle
+        
+        
+     artio_file artio_fileset_open(char * file_prefix, int type, artio_context context) 
+
+     void wrap_artio_fileset_open(char *file_prefix, int type)
+   
+
+#python only passes numbers and strings... no pointers or structures
+cpdef read_header(char * file_prefix, int type): 
+     wrap_artio_fileset_open(file_prefix, type)
+        
+

yt/frontends/artio/api.py

 """
-API for yt.frontends.ramses
+API for yt.frontends.artio
 
 Author: Matthew Turk <matthewturk@gmail.com>
 Affiliation: UCSD
 """
 
 from .data_structures import \
-      RAMSESStaticOutput
+      ARTIOStaticOutput
 
 from .fields import \
-      RAMSESFieldInfo, \
-      add_ramses_field
+      ARTIOFieldInfo, \
+      add_artio_field
 
 from .io import \
-      IOHandlerRAMSES
+      IOHandlerARTIO

yt/frontends/artio/artio_caller.pyx

-#struct artio_context_struct {
-#        MPI_Comm comm;
-#};
-cdef struct artio_context_struct: #in artio_mpi.h MPI_Comm
-    int comm
-
-#typedef struct artio_context_struct * artio_context;
-ctypedef artio_context_struct * artio_context
-
-cdef extern from "artio_headers/artio.c":
-
-     ctypedef int int64_t 
-
-     ctypedef struct param: 
-         int key_length
-         char key[64]
-         int val_length
-         int type 
-         char *value
-         
-     ctypedef struct list: 
-         param * head
-         param * tail
-         param * cursor
-         int iterate_flag 
-         int endian_swap 
-
-     #    artio_internal.h
-     ctypedef struct artio_particle_file:
-        #	artio_fh *ffh
-        int num_particle_files
-        int64_t *file_sfc_index
-        int64_t cache_sfc_begin
-        int64_t cache_sfc_end
-        int64_t *sfc_offset_table
-        
-        #/* maintained for consistency and user-error detection */
-        int num_species
-        int cur_file
-        int cur_species
-        int cur_particle
-        int64_t cur_sfc
-        int *num_primary_variables
-        int *num_secondary_variables
-        int *num_particles_per_species
-        
-     ctypedef struct artio_grid_file:
-        #        artio_fh *ffh
-        int num_grid_variables
-        int num_grid_files
-        int64_t *file_sfc_index
-        int64_t cache_sfc_begin
-        int64_t cache_sfc_end
-        int64_t *sfc_offset_table
-        
-        int file_max_level
-        #/* maintained for consistency and user-error detection */
-        int cur_file
-        int cur_num_levels
-        int cur_level
-        int cur_octs
-        int64_t cur_sfc
-        int *octs_per_level
-        
-        
-     ctypedef struct artio_file:
-        char file_prefix[256]
-        int endian_swap
-        int open_type
-        int open_mode
-        int rank
-        int num_procs
-        artio_context context #artio_mpi.h MPI_Comm
-        
-        int64_t *proc_sfc_index
-        int64_t proc_sfc_begin
-        int64_t proc_sfc_end
-        int64_t num_root_cells
-        
-        list param_list
-        artio_grid_file grid
-        artio_particle_file particle
-        
-        
-     artio_file artio_fileset_open(char * file_prefix, int type, artio_context context) 
-
-     void wrap_artio_fileset_open(char *file_prefix, int type)
-   
-
-#python only passes numbers and strings... no pointers or structures
-cpdef read_header(char * file_prefix, int type): 
-     wrap_artio_fileset_open(file_prefix, type)
-        
-

yt/frontends/artio/artio_headers/ARTIO_amr_data.hh

+/*
+		This file is part of libARTIO++ 
+			a C++ library to access snapshot files 
+			generated by the simulation code ARTIO by R. Teyssier
+		
+    Copyright (C) 2008-09  Oliver Hahn, ojha@gmx.de
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __ARTIO_AMR_DATA_HH
+#define __ARTIO_AMR_DATA_HH
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <iterator>
+
+#include "FortranUnformatted_IO.hh"
+#include "ARTIO_info.hh"
+#include "ARTIO_amr_data.hh"
+
+#ifndef FIX
+#define FIX(x)	((int)((x)+0.5))
+#endif
+
+#define ACC_NL(cpu,lvl)   ((cpu)+m_header.ncpu*(lvl))
+#define ACC_NLIT(cpu,lvl) ((cpu)+m_plevel->m_header.ncpu*(lvl))
+#define ACC_NB(cpu,lvl)   ((cpu)+m_header.nboundary*(lvl))
+
+#define LENGTH_POINTERLISTS 4096
+#define ENDPOINT ((unsigned)(-1))
+
+namespace ARTIO{
+namespace AMR{
+	
+/**************************************************************************************\
+ *** auxiliary datatypes **************************************************************
+\**************************************************************************************/
+	
+template< typename real_t>
+struct vec{
+	real_t x,y,z;
+		
+	vec( real_t x_, real_t y_, real_t z_ )
+	: x(x_),y(y_),z(z_)
+	{ }
+		
+	vec( const vec& v )
+	: x(v.x),y(v.y),z(v.z)
+	{ }
+		
+	vec( void )
+	: x(0.0), y(0.0), z(0.0)
+	{ }
+};
+
+/**************************************************************************************\
+ *** AMR cell base types **************************************************************
+\**************************************************************************************/
+
+template <typename id_t=unsigned, typename real_t=float>
+class cell_locally_essential{
+public:
+	id_t m_neighbour[6];
+	id_t m_father;
+	id_t m_son[8];
+	real_t m_xg[3];
+	id_t m_cpu;
+
+	char m_pos;
+	
+	cell_locally_essential(){}
+		
+	bool is_refined( int ison ) const
+	{	return ((int)m_son[ison]!=-1); }	
+};
+
+template <typename id_t=unsigned, typename real_t=float>
+class cell_simple{
+public:
+	id_t m_son[1];
+	real_t m_xg[3];
+	id_t m_cpu;
+
+	char m_pos;
+	
+	cell_simple(){}
+		
+	bool is_refined( int ison=0 ) const
+	{	return ((int)m_son[0]!=-1); }	
+};
+
+
+//.... some type traits that are used to distinguish what data needs to be read ....//
+template<class X> struct is_locally_essential
+{ enum { check=false }; };
+
+template<> struct is_locally_essential<cell_locally_essential<> > 
+{ enum { check=true }; };
+
+template<> struct is_locally_essential<cell_simple<> > 
+{ enum { check=false }; };
+
+
+/**************************************************************************************\
+ *** AMR level class definition, subsumes a collection of AMR cells *******************
+\**************************************************************************************/
+
+//! AMR level implementation
+template< typename Cell_ >
+class level{
+public:
+	unsigned m_ilevel;
+	std::vector< Cell_ > m_level_cells;
+	
+	double m_xc[8];			//!< relative x-offsets of the 8 children for current refinement level
+	double m_yc[8];			//!< relative y-offsets of the 8 children for current refinement level
+	double m_zc[8];			//!< relative z-offsets of the 8 children for current refinement level
+	
+	typedef typename std::vector< Cell_ >::iterator iterator;
+	typedef typename std::vector< Cell_ >::const_iterator const_iterator;
+
+	double 		m_dx;
+	unsigned 	m_nx;
+
+	level( unsigned ilevel )
+	: m_ilevel( ilevel )
+	{
+		m_dx = pow(0.5,ilevel+1);
+		m_nx = (unsigned)(1.0/m_dx);
+		
+		for( unsigned k=1; k<=8; k++ )
+		{
+			//... initialize positions of son cell centres
+			//... relative to parent cell centre
+			unsigned 
+				iz=(k-1)/4,
+				iy=(k-1-4*iz)/2,
+				ix=(k-1-2*iy-4*iz);
+			m_xc[k-1]=((double)ix-0.5f)*m_dx;
+			m_yc[k-1]=((double)iy-0.5f)*m_dx;
+			m_zc[k-1]=((double)iz-0.5f)*m_dx;
+		}
+	}
+	
+	void register_cell( const Cell_ & cell )
+	{ m_level_cells.push_back( cell ); }
+	
+	const_iterator begin( void ) const{ return m_level_cells.begin(); }
+	iterator begin( void ){ return m_level_cells.begin(); }
+	
+	const_iterator end( void ) const{ return m_level_cells.end(); }
+	iterator end( void ){ return m_level_cells.end(); }
+	
+	Cell_& operator[]( unsigned i )
+	{ return m_level_cells[i]; }
+	
+	unsigned size( void ) { return m_level_cells.size(); }
+};
+
+
+/**************************************************************************************\
+ *** constants ************************************************************************
+\**************************************************************************************/
+
+//! neighbour cell access pattern
+const static int nbor_cell_map[6][8] = 
+{ 
+	{1,0,3,2,5,4,7,6},
+	{1,0,3,2,5,4,7,6},
+	{2,3,0,1,6,7,4,5},
+	{2,3,0,1,6,7,4,5},
+	{4,5,6,7,0,1,2,3},
+	{4,5,6,7,0,1,2,3} 
+};
+
+
+/**************************************************************************************\
+ *** AMR tree class implements the hierarchy of AMR levels with links *****************
+\**************************************************************************************/
+
+/*!
+ * @class tree
+ * @brief encapsulates the hierarchical AMR structure data from a ARTIO simulation snapshot
+ *
+ * This class provides low-level read access to ARTIO amr_XXXXX.out files. 
+ * Data from a given list of computational domains can be read and is
+ * stored in internal datastructures.
+ * Access to hydrodynamical variables stored in the cells is provided 
+ * through member functions of class ARTIO_hydro_level and iterators 
+ * provided by this class
+ * @sa ARTIO_hydro_level
+ */
+template< class Cell_, class Level_ >
+class tree
+{
+
+public:
+	
+	//! header amr meta-data structure, for details see also ARTIO source code (file amr/init_amr.f90)
+	struct header{ 
+		std::vector< int > nx;			//!< base mesh resolution [3-vector]
+		std::vector< int > nout;		//!< 3 element vector: [noutput2, iout2, ifout2]
+		std::vector< int > nsteps;		//!< 2 element vector: [nstep, nstep_coarse]
+		int ncpu;						//!< number of CPUs (=computational domains) in the simulation
+		int ndim;						//!< number of spatial dimensions
+		int nlevelmax;					//!< maximum refinement level allowed
+		int ngridmax;					//!< maxmium number of grid cells stored per CPU
+		int nboundary;					//!< number of boundary cells (=ghost cells)
+		int ngrid_current;				//!< currently active grid cells
+		double boxlen;					//!< length of the simulation box
+		std::vector<double> tout;		//!< output times (1..noutput)
+		std::vector<double> aout;		//!< output times given as expansion factors (1..noutput)
+		std::vector<double> dtold;		//!< old time steps (1..nlevelmax)
+		std::vector<double> dtnew;		//!< next time steps (1..nlevelmax)
+		std::vector<double> stat;		//!< some diagnostic snapshot meta data: [const, mass_tot0, rho_tot]
+		std::vector<double> cosm;		//!< cosmological meta data: [omega_m, omega_l, omega_k, omega_b, h0, aexp_ini, boxlen_ini)
+		std::vector<double> timing;		//!< timing information: [aexp, hexp, aexp_old, epot_tot_int, epot_tot_old]
+		double t;						//!< time stamp of snapshot
+		double mass_sph;				//!< mass threshold used to flag for cell refinement
+	};
+	
+	std::vector< Level_ > m_AMR_levels;	//! STL vector holding the AMR level hierarchy
+	
+	std::vector<unsigned> 
+		m_headl,						//!< head indices, point to first active cell
+		m_numbl, 						//!< number of active cells
+		m_taill;						//!< tail indices, point to last active cell
+
+	int m_cpu;							//! index of computational domain being accessed
+	int m_minlevel;						//! lowest refinement level to be loaded
+	int m_maxlevel;						//! highest refinement level to be loaded
+	std::string m_fname;				//! the snapshot filename amr_XXXXX.out
+	unsigned m_ncoarse;					//! number of coarse grids
+	struct header m_header;				//! the header meta data
+	
+	
+protected:
+	
+  //! prototypical grid iterator
+  /*! iterates through cells on one level, provides access to neighbours, parents, children
+   */
+	template< typename TreePointer_, typename Index_=unsigned >
+	class proto_iterator{
+		public:
+		
+			friend class tree;
+			
+			typedef Index_ value_type;
+			typedef Index_& reference;
+			typedef Index_* pointer;
+			
+		protected:
+		
+			Index_
+				m_ilevel,		//!< refinement level on which we iteratre
+				m_icpu;			//!< domain on which we iterate
+				
+			typedef typename std::vector< Cell_ >::const_iterator Iterator_;
+			Iterator_	 m_cell_it; //!< basis iterator that steps through cells on one level
+			TreePointer_ m_ptree; //!< pointer to associated tree object
+			
+			
+			//! low_level construtor that should only be used from within AMR_tree
+			proto_iterator( unsigned level_, unsigned cpu_, Iterator_ it_, TreePointer_ ptree_ )
+			: m_ilevel(level_), m_icpu(cpu_), m_cell_it(it_), m_ptree( ptree_ )
+			{ }
+			
+		public:
+		
+			//! this is either the copy-constructor or a constructor for implicit type conversion
+			template< typename X >
+			proto_iterator( proto_iterator<X> &it )
+			: m_ilevel( it.m_ilevel ), m_icpu( it.m_icpu ), m_cell_it( it.m_cell_it ), m_ptree( it.m_ptree )
+			{ }
+	
+			//! empty constructor, doesn't initialize anything
+			proto_iterator()
+			: m_ilevel(0), m_icpu(0), m_ptree(NULL)
+			{ }
+			
+			//! test for equality between two proto_iterator instantiations
+			template< typename X >
+			inline bool operator==( const proto_iterator<X>& x ) const
+			{ return m_cell_it==x.m_cell_it; }
+
+			//! test for inequality between two proto_iterator instantiations	
+			template< typename X >
+			inline bool operator!=( const proto_iterator<X>& x ) const
+			{ return m_cell_it!=x.m_cell_it; }
+
+			//! iterate forward, prefix
+			inline proto_iterator& operator++()
+			{ ++m_cell_it; return *this; }
+			
+			//! iterate forward, postfix
+			inline proto_iterator operator++(int)
+			{ proto_iterator<TreePointer_> it(*this); operator++(); return it; }
+
+            inline void next(void) { operator++(); }
+			
+			//! iterate backward, prefix
+			inline proto_iterator& operator--()
+			{ --m_cell_it; return *this; }
+			
+			//! iterate backward, postfix
+			inline proto_iterator operator--(int)
+			{ proto_iterator<TreePointer_> it(*this); operator--(); return it; }
+			
+			//! iterate several forward
+			inline proto_iterator operator+=(int i)
+			{ proto_iterator<TreePointer_> it(*this); m_cell_it+=i; return it; }
+			
+			//! iterate several backward
+			inline proto_iterator operator-=(int i)
+			{ proto_iterator<TreePointer_> it(*this); m_cell_it-=i; return it; }
+			
+			//! assign two proto_iterators, this will fail if no typecast between X and TreePoint_ exists
+			template< typename X >
+			inline proto_iterator& operator=(const proto_iterator<X>& x)
+			{ m_cell_it = x.m_cell_it; m_ilevel = x.m_ilevel; m_icpu = x.m_icpu; m_ptree = x.m_ptree; return *this; }
+			
+			//! iterator dereferencing, returns an array index
+			inline Cell_ operator*() const
+			{ return *m_cell_it; }
+			
+            inline Index_ get_cell_father() const { return (*m_cell_it).m_father; }
+            inline bool is_finest(int ison) { return !((*m_cell_it).is_refined(ison)); }
+			
+			//! move iterator to a child grid
+			/*!
+			 * @param  ind index of the child grid in the parent oct (0..7)
+			 * @return iterator pointing to child grid if it exists, otherwise 'end' of currrent level
+			 */
+			//template< typename X >
+			inline proto_iterator& to_child( int ind )
+			{
+				if( !(*m_cell_it).is_refined(ind) )
+					 return (*this = m_ptree->end( m_ilevel ));
+				++m_ilevel;
+				m_cell_it = m_ptree->m_AMR_levels[m_ilevel].begin()+(*m_cell_it).m_son[ind];
+				return *this;
+			}
+			
+			
+			//! get an iterator to a child grid
+			/*!
+			 * @param  ind index of the child grid in the parent oct (0..7)
+			 * @return iterator pointing to child grid if it exists, otherwise 'end' of currrent level
+			 */
+			//template< typename X >
+			inline proto_iterator get_child( int ind ) const
+			{
+				proto_iterator it(*this);
+				it.to_child( ind );
+				return it;
+			}
+			
+			inline char get_child_from_pos( double x, double y, double z ) const
+			{
+				bool  bx,by,bz;
+				bx = x > (*m_cell_it).m_xg[0];
+				by = y > (*m_cell_it).m_xg[1];
+				bz = z > (*m_cell_it).m_xg[2];
+				
+				//std::cerr << "(" << bx << ", " << by << ", " << bz << ")\n";
+				
+				return (char)bx+2*((char)by+2*(char)bz);
+			}
+			
+			
+			//! move iterator to the parent grid
+			/*! 
+			 * @return iterator pointing to the parent grid if it exists, 'end' of the current level otherwise
+			 */
+			
+			inline proto_iterator& to_parent( void )
+			{
+				if( m_ilevel==0 )
+					return (*this = m_ptree->end( m_ilevel ));
+				--m_ilevel;
+				m_cell_it = m_ptree->m_AMR_levels[m_ilevel].begin()+(*m_cell_it).m_father;
+				return *this;
+			}
+			
+			//! query whether a given child cell is refined
+			inline bool is_refined( int i ) const
+			{
+				return (*m_cell_it).is_refined(i);
+			}
+			
+			//! get an iterator to the parent grid
+			/*! 
+			 * @return iterator pointing to the parent grid if it exists, 'end' of the current level otherwise
+			 */
+			inline proto_iterator get_parent( void ) const
+			{
+				proto_iterator it(*this);
+				it.to_parent();
+				return it;
+			}
+			
+			//! move iterator to spatial neighbour grid
+			/*!
+			 * @param ind index of neighbour (0..5)
+			 * @return iterator pointing to neighbour grid if it exists, otherwise 'end' of currrent level
+			 */
+			inline proto_iterator& to_neighbour( int ind )
+			{
+				unsigned icell = nbor_cell_map[ind][(int)(*m_cell_it).m_pos];
+				m_cell_it = m_ptree->m_AMR_levels[m_ilevel-1].begin()+(*m_cell_it).m_neighbour[ind];
+				
+				if( !(*m_cell_it).is_refined(icell) )
+					return (*this = m_ptree->end(m_ilevel));
+				
+				m_cell_it = m_ptree->m_AMR_levels[m_ilevel].begin()+(*m_cell_it).m_son[icell];
+				return *this;
+			}
+			
+			//! get an iterator to spatial neighbour grid
+			/*!
+			 * @param ind index of neighbour (0..5)
+			 * @return iterator pointing to neighbour grid if it exists, otherwise 'end' of currrent level
+			 */
+			inline proto_iterator& get_neighbour( int ind )
+			{
+				proto_iterator it(*this);
+				it.to_neighbour(ind);
+				return it;
+			}
+			
+			inline Index_ get_level( void ) const
+			{ return m_ilevel; }
+			
+			inline int get_domain( void ) const
+			{ return (*m_cell_it).m_cpu; }
+			
+			inline int get_absolute_position( void ) const 
+			{
+				return (unsigned)(std::distance<Iterator_>(m_ptree->m_AMR_levels[m_ilevel].begin(),m_cell_it));
+			}
+			
+			
+	};
+	
+public:
+	
+	typedef proto_iterator<const tree*> const_iterator;
+	typedef proto_iterator<tree*>       iterator;
+	
+
+
+protected:
+	
+	//! read header meta data from amr snapshot file
+	void read_header( void );
+	
+
+	
+	//! generate the amr_XXXXX.out filename for a given computational domain
+	/*! @param icpu index of comutational domain (base 1)
+	 */
+	std::string gen_fname( int icpu );
+	
+	//! generate the amr_XXXXX.out filename from the path to the info_XXXXX.out file
+	std::string rename_info2amr( const std::string& info );
+
+
+	#define R_SQR(x) ((x)*(x))
+	
+	template< typename Real_ >
+	inline bool ball_intersection( const vec<Real_>& xg, double dx2, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dmin = 0, bmin, bmax;
+		
+		//.. x ..//
+		bmin = xg.x-dx2;
+		bmax = xg.x+dx2;
+		if( xc.x < bmin ) dmin += R_SQR(xc.x - bmin ); else
+		if( xc.x > bmax ) dmin += R_SQR(xc.x - bmax );
+			
+		//.. y ..//
+		bmin = xg.y-dx2;
+		bmax = xg.y+dx2;
+		if( xc.y < bmin ) dmin += R_SQR(xc.y - bmin ); else
+		if( xc.y > bmax ) dmin += R_SQR(xc.y - bmax );
+		
+		//.. x ..//
+		bmin = xg.z-dx2;
+		bmax = xg.z+dx2;
+		if( xc.z < bmin ) dmin += R_SQR(xc.z - bmin ); else
+		if( xc.z > bmax ) dmin += R_SQR(xc.z - bmax );
+		
+		if( dmin <= r2 ) return true;
+		return false;
+	}
+	
+	template< typename Real_ >
+	inline bool shell_intersection( const vec<Real_>& xg, double dx2, const vec<Real_>& xc, Real_ r1_2, Real_ r2_2 )
+	{
+		Real_ dmax = 0, dmin = 0, a, b, bmin, bmax;
+		if( r1_2 > r2_2 ) std::swap(r1_2,r2_2);
+			
+		//.. x ..//
+		bmin = xg.x-dx2;
+		bmax = xg.x+dx2;
+		a = R_SQR( xc.x - bmin );
+		b = R_SQR( xc.x - bmax );
+		dmax += std::max( a, b );
+		if( xc.x < bmin ) dmin += a; else
+		if( xc.x > bmax ) dmin += b;
+		
+		//.. y ..//
+		bmin = xg.y-dx2;
+		bmax = xg.y+dx2;
+		a = R_SQR( xc.y - bmin );
+		b = R_SQR( xc.y - bmax );
+		dmax += std::max( a, b );
+		if( xc.y < bmin ) dmin += a; else
+		if( xc.y > bmax ) dmin += b;
+			
+		//.. z ..//
+		bmin = xg.z-dx2;
+		bmax = xg.z+dx2;
+		a = R_SQR( xc.z - bmin );
+		b = R_SQR( xc.z - bmax );
+		dmax += std::max( a, b );
+		if( xc.z < bmin ) dmin += a; else
+		if( xc.z > bmax ) dmin += b;
+		
+		
+		if( dmin <= r2_2 && r1_2 <= dmax ) return true;
+		return false;
+	}
+	
+	template< typename Real_ >
+	inline bool sphere_intersection( const vec<Real_>& xg, double dx2, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dmax = 0, dmin = 0, a, b, bmin, bmax;
+		
+		//.. x ..//
+		bmin = xg.x-dx2;
+		bmax = xg.x+dx2;
+		a = R_SQR( xc.x - bmin );
+		b = R_SQR( xc.x - bmax );
+		dmax += std::max( a, b );
+		if( xc.x < bmin ) dmin += a; else
+		if( xc.x > bmax ) dmin += b;
+		
+		//.. y ..//
+		bmin = xg.y-dx2;
+		bmax = xg.y+dx2;
+		a = R_SQR( xc.y - bmin );
+		b = R_SQR( xc.y - bmax );
+		dmax += std::max( a, b );
+		if( xc.y < bmin ) dmin += a; else
+		if( xc.y > bmax ) dmin += b;
+			
+		//.. z ..//
+		bmin = xg.z-dx2;
+		bmax = xg.z+dx2;
+		a = R_SQR( xc.z - bmin );
+		b = R_SQR( xc.z - bmax );
+		dmax += std::max( a, b );
+		if( xc.z < bmin ) dmin += a; else
+		if( xc.z > bmax ) dmin += b;
+		
+		
+		if( dmin <= r2 && r2 <= dmax ) return true;
+		return false;
+	}
+	
+	#undef R_SQR
+	
+public:
+	
+	//! low-level constructor - should not be called from outside because then you can really screw up things
+	/*!
+	 * @param snap the associated ARTIO::snapshot object
+	 * @param cpu domain for which to read the AMR tree
+	 * @param maxlevel maximum refinement level to consider
+	 * @param minlevel minimum refinement level to consider (default=1)
+	 */
+	tree( ARTIO::snapshot& snap, int cpu, int maxlevel, int minlevel=1 )
+	: m_cpu( cpu ), m_minlevel( minlevel ), m_maxlevel( maxlevel ), m_fname( rename_info2amr(snap.m_filename) )
+	{ 
+		read_header();
+
+		if( cpu > m_header.ncpu || cpu <= 0)
+			throw std::runtime_error("ARTIO_particle_data: expect to read from out of range CPU.");
+		
+	}
+	
+	//! perform the read operation of AMR data
+	void read( void );
+	
+	//! end const_iterator for given refinement level
+	const_iterator end( int ilevel ) const
+	{ 
+		if( ilevel <= m_maxlevel )
+			return const_iterator( ilevel, m_cpu, m_AMR_levels.at(ilevel).end(), this );
+		else
+			return const_iterator( ilevel, m_cpu, m_AMR_levels.at(0).end(), this );
+	}
+	
+	//! end iterator for given refinement level
+	iterator end( int ilevel )
+	{ 
+		if( ilevel <= m_maxlevel )
+			return iterator( ilevel, m_cpu, m_AMR_levels.at(ilevel).end(), this ); 
+		else
+			return iterator( ilevel, m_cpu, m_AMR_levels.at(0).end(), this ); 
+	}
+	
+	//! begin const_iterator for given refinement level
+	const_iterator begin( int ilevel ) const
+	{	
+		if( ilevel <= m_maxlevel )
+			return const_iterator( ilevel, m_cpu, m_AMR_levels.at(ilevel).begin(), this ); 
+		else
+			return this->end(ilevel);
+	}
+	
+	//! begin iterator for given refinement level
+	iterator begin( int ilevel )
+	{	
+		if( ilevel <= m_maxlevel )
+			return iterator( ilevel, m_cpu, m_AMR_levels.at(ilevel).begin(), this ); 
+		else
+			return this->end(ilevel);
+	}
+	
+	
+	//! return the center of a child cell associated with a grid iterator
+	/*!
+	 * @param it grid iterator
+	 * @param ind sub-cell index (0..7)
+	 * @return vec vector containing the coordinates
+	 */
+	template< typename Real_ >
+	inline vec<Real_> cell_pos( const iterator& it, unsigned ind )
+	{
+		vec<Real_> pos;
+		pos.x = (*it).m_xg[0]+m_AMR_levels[it.m_ilevel].m_xc[ind];
+		pos.y = (*it).m_xg[1]+m_AMR_levels[it.m_ilevel].m_yc[ind];
+		pos.z = (*it).m_xg[2]+m_AMR_levels[it.m_ilevel].m_zc[ind];
+		return pos;
+	}
+	
+	//! return the center of the grid associated with a grid iterator
+	/*!
+	 * @param it grid iterator
+	 * @return vec vector containing the coordinates
+	 */
+	template< typename Real_ >
+	inline vec<Real_> grid_pos( const iterator& it )
+	{
+		vec<Real_> pos;
+		pos.x = (*it).m_xg[0];
+		pos.y = (*it).m_xg[1];
+		pos.z = (*it).m_xg[2];
+		return pos;
+	}
+	
+	template< typename Real_ >
+	inline bool ball_intersects_grid( const iterator& it, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level());
+		vec<Real_> xg = grid_pos<Real_>(it);
+		return ball_intersection( xg, dx2, xc, r2 );
+	}
+	
+	template< typename Real_ >
+	inline bool ball_intersects_cell( const iterator& it, char ind, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level()+1);
+		vec<Real_> xg = cell_pos<Real_>(it,ind);
+		return ball_intersection( xg, dx2, xc, r2 );
+	}
+		
+	template< typename Real_ >
+	inline bool shell_intersects_grid( iterator& it, const vec<Real_>& xc, Real_ r1_2, Real_ r2_2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level());
+		vec<Real_> xg = grid_pos<Real_>(it);
+		return shell_intersection( xg, dx2, xc, r1_2, r2_2 );
+	}
+	
+	template< typename Real_ >
+	inline bool shell_intersects_cell( iterator& it, char ind, const vec<Real_>& xc, Real_ r1_2, Real_ r2_2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level()+1);
+		vec<Real_> xg = cell_pos<Real_>(it,ind);
+		return shell_intersection( xg, dx2, xc, r1_2, r2_2 );
+	}
+		
+	template< typename Real_ >
+	inline bool sphere_intersects_grid( const iterator& it, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level());
+		vec<Real_> xg = grid_pos<Real_>(it);
+		return sphere_intersection( xg, dx2, xc, r2 );
+	}
+	
+	template< typename Real_ >
+	inline bool sphere_intersects_cell( const iterator& it, char ind, const vec<Real_>& xc, Real_ r2 )
+	{
+		Real_ dx2 = 0.5/pow(2,it.get_level()+1);
+		vec<Real_> xg = cell_pos<Real_>(it,ind);
+		return sphere_intersection( xg, dx2, xc, r2 );
+	}
+	
+};
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+template< class Cell_, class Level_ >
+void tree<Cell_,Level_>::read_header( void )
+{
+    FortranUnformatted ff( gen_fname(m_cpu) );
+    std::cout << "snl: HI I AM HERE AND ALIVE!!!!!!";
+	
+	//-- read header data --//
+	
+	ff.read( m_header.ncpu );
+	ff.read( m_header.ndim );
+	ff.read<unsigned>( std::back_inserter(m_header.nx) );
+	ff.read( m_header.nlevelmax );
+	ff.read( m_header.ngridmax );
+	ff.read( m_header.nboundary );
+	ff.read( m_header.ngrid_current );
+	ff.read( m_header.boxlen );
+	
+	ff.read<unsigned>( std::back_inserter(m_header.nout) );
+	ff.read<double>( std::back_inserter(m_header.tout) );
+	ff.read<double>( std::back_inserter(m_header.aout) );
+	ff.read( m_header.t );
+	ff.read<double>( std::back_inserter(m_header.dtold) );
+	ff.read<double>( std::back_inserter(m_header.dtnew) );
+	ff.read<unsigned>( std::back_inserter(m_header.nsteps) );
+	ff.read<double>( std::back_inserter(m_header.stat) );
+	ff.read<double>( std::back_inserter(m_header.cosm) );
+	ff.read<double>( std::back_inserter(m_header.timing) );
+	ff.read( m_header.mass_sph );
+	
+	m_ncoarse = m_header.nx[0]*m_header.nx[1]*m_header.nx[2];
+}
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+template< class Cell_, class Level_ >
+std::string tree<Cell_,Level_>::gen_fname( int icpu )
+{
+	std::string fname;
+	char ext[32];
+	fname = m_fname;
+	fname.erase(fname.rfind('.')+1);
+	sprintf(ext,"out%05d",icpu);
+	fname.append(std::string(ext));
+	return fname;
+}
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+template< class Cell_, class Level_ >
+std::string tree<Cell_,Level_>::rename_info2amr( const std::string& info )
+{
+	std::string amr;
+	unsigned ii = info.rfind("info");
+	amr = info.substr(0,ii)+"amr" + info.substr(ii+4, 6) + ".out00001";
+	return amr;
+}
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+template< class Cell_, class Level_ >
+void tree<Cell_,Level_>::read( void )
+{
+	// indexing map used to associate ARTIO internal cell IDs with new IDs
+	std::map<unsigned,unsigned> m_ind_grid_map;
+
+	std::vector<int> cell_cpu;
+	std::vector<unsigned> cell_level;
+	std::vector<unsigned> itmp;
+	
+	
+	typename std::vector< Level_ >::iterator amr_level_it;
+	
+	FortranUnformatted ff( gen_fname( m_cpu ) );
+		
+	//.. skip header entries ..//
+	ff.skip_n_from_start( 19 ); 			//.. skip header 
+		
+	//+ headl + taill
+	ff.read<int>( std::back_inserter(m_headl) );
+	ff.read<int>( std::back_inserter(m_taill) );
+	ff.read<int>( std::back_inserter(m_numbl) );
+	
+	//.. skip numbtot
+	ff.skip_n( 1 ); 						
+		
+	std::vector<int> ngridbound;
+	if( m_header.nboundary > 0 ){
+		ff.skip_n( 2 ); 					//.. skip headb and tailb 
+		ff.read<int>( std::back_inserter(ngridbound) ); 				//.. read numbb
+	}
+		
+	ff.skip_n( 6 ); //..skip (1)free_mem+(2)ordering+(3)bound_key+
+					//..     (4)coarse_son+(5)coarse_flag1+(6)coarse_cpu_map
+	
+	
+	if( /*m_minlevel < 1 ||*/ m_maxlevel > m_header.nlevelmax || m_minlevel > m_maxlevel )
+		throw std::runtime_error("ARTIO_amr_level::read_level : requested level is invalid.");
+		
+
+	m_ind_grid_map.insert( std::pair<unsigned,unsigned>(0,ENDPOINT) );
+	
+	FortranUnformatted::streampos spos = ff.tellg();
+	
+	m_minlevel = 0;
+	
+	//... create indexing map ...//
+	for( int ilvl = 0; ilvl<=std::min(m_maxlevel+1, m_header.nlevelmax); ++ilvl ){
+		unsigned gridoff = 0;
+		for( int icpu=0; icpu<m_header.ncpu+m_header.nboundary; ++icpu ){
+			if( icpu < m_header.ncpu && m_numbl[ACC_NL(icpu,ilvl)] == 0 )
+				continue;
+			else if( icpu >= m_header.ncpu && ngridbound[ACC_NB(icpu-m_header.ncpu,ilvl)] == 0 )
+				continue;
+				
+			if( ilvl >= m_minlevel ){
+				std::vector<int> ind_grid;
+				ff.read<int>( std::back_inserter(ind_grid) );
+				for( unsigned i=0; i<ind_grid.size(); ++i ){
+					m_ind_grid_map.insert( std::pair<unsigned,unsigned>(ind_grid[i],gridoff++) );
+				}
+				ind_grid.clear();
+			}else
+				ff.skip();
+				
+			ff.skip_n( 3+3+6+8+8+8 );
+		}
+		if( ff.eof() ){
+			//std::cerr << "eof reached in fortran read operation\n";
+			m_maxlevel = ilvl;//+1;
+			break;
+		}
+	}
+	
+	ff.seekg( spos );
+	
+	m_AMR_levels.clear();
+	
+	 
+	//... loop over levels ...//
+	for( int ilvl = 0; ilvl<=m_maxlevel; ++ilvl ){
+		m_AMR_levels.push_back( Level_(ilvl) );
+		Level_ &currlvl = m_AMR_levels.back();
+			
+		for( int icpu=0; icpu<m_header.ncpu+m_header.nboundary; ++icpu ){
+			if( icpu < m_header.ncpu && m_numbl[ACC_NL(icpu,ilvl)] == 0 )
+				continue;
+			else if( icpu >= m_header.ncpu && ngridbound[ACC_NB(icpu-m_header.ncpu,ilvl)] == 0 )
+				continue;
+			
+			if( ilvl >= m_minlevel ){
+				unsigned gridoff = currlvl.size();
+			
+				std::vector<int> ind_grid;
+				ff.read<int>( std::back_inserter(ind_grid) );
+				for( unsigned i=0; i<ind_grid.size(); ++i ){
+					currlvl.register_cell( Cell_() );
+					//.. also set owning cpu in this loop...
+					currlvl[ i+gridoff ].m_cpu = icpu+1;
+				}
+			
+				//... pointers to next and previous octs ..//
+				ff.skip();
+				ff.skip();
+			
+				//... oct x-coordinates ..//
+				std::vector<float> ftmp;
+				ff.read<double>( std::back_inserter(ftmp) );
+				for( unsigned j=0; j<ftmp.size(); ++j )
+					currlvl[ j+gridoff ].m_xg[0] = ftmp[j];
+				ftmp.clear();
+			
+				//... oct y-coordinates ..//
+				ff.read<double>( std::back_inserter(ftmp) );
+				for( unsigned j=0; j<ftmp.size(); ++j )
+					currlvl[ j+gridoff ].m_xg[1] = ftmp[j];
+				ftmp.clear();
+			
+				//... oct y-coordinates ..//
+				ff.read<double>( std::back_inserter(ftmp) );
+				for( unsigned j=0; j<ftmp.size(); ++j )
+					currlvl[ j+gridoff ].m_xg[2] = ftmp[j];
+				ftmp.clear();
+			
+			
+				//... father indices
+				if( is_locally_essential<Cell_>::check ){
+					ff.read<int>( std::back_inserter(itmp) ); 
+					for( unsigned j=0; j<itmp.size(); ++j ){
+						currlvl[ j+gridoff ].m_pos    = (itmp[j]-m_ncoarse-1)/m_header.ngridmax;
+						currlvl[ j+gridoff ].m_father = m_ind_grid_map[ (itmp[j]-m_ncoarse)%m_header.ngridmax ];
+					}
+					itmp.clear();
+				}else
+					ff.skip();
+				
+				
+				//... neighbour grids indices
+				if( is_locally_essential<Cell_>::check )
+				{
+					for( unsigned k=0; k<6; ++k ){
+						ff.read<int>( std::back_inserter(itmp) );
+						for( unsigned j=0; j<itmp.size(); ++j )
+							currlvl[j+gridoff].m_neighbour[k] = m_ind_grid_map[ (itmp[j]-m_ncoarse)%m_header.ngridmax ];
+						itmp.clear();
+					}
+				}else
+					ff.skip_n( 6 );
+					
+			
+				//... son cell indices
+				for( unsigned ind=0; ind<8; ++ind ){
+					ff.read<int>( std::back_inserter(itmp) );
+					for( unsigned k=0; k<itmp.size(); ++k ){
+						currlvl[k+gridoff].m_son[ind] = m_ind_grid_map[ itmp[k] ];
+					}
+					itmp.clear();
+				}
+
+				//.. skip cpu + refinement map
+				ff.skip_n( 8+8 );
+			}else{
+				//...skip entire record
+				ff.skip_n( 3+3+1+6+8+8+8 );		
+			}
+		}
+	}
+}
+
+} //namespace AMR
+} //namespace ARTIO
+
+
+
+#undef FIX
+
+#endif //__ARTIO_AMR_DATA_HH

yt/frontends/artio/artio_headers/ARTIO_hydro_data.hh

+/*
+		This file is part of libARTIO++
+			a C++ library to access snapshot files
+			generated by the simulation code ARTIO by R. Teyssier
+
+    Copyright (C) 2008-09  Oliver Hahn, ojha@gmx.de
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __ARTIO_HYDRO_DATA_HH
+#define __ARTIO_HYDRO_DATA_HH
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+
+#include "FortranUnformatted_IO.hh"
+#include "ARTIO_info.hh"
+#include "ARTIO_amr_data.hh"
+
+namespace ARTIO{
+namespace HYDRO{
+
+//! internal hydro variable indices
+enum hydro_var
+{
+	density     = 1,
+	velocity_x  = 2,
+	velocity_y  = 3,
+	velocity_z  = 4,
+	pressure    = 5,
+	metallicity = 6
+};
+
+//! names of possible variables stored in a ARTIO hydro file
+const char artio_hydro_variables[][64] = {
+	{"Density"},
+	{"x-velocity"},
+	{"y-velocity"},
+	{"z-velocity"},
+	{"Pressure"},
+	{"Metallicity"} };
+
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+/*!
+ * @class ARTIO::HYDRO::proto_data
+ * @brief base class for all things hydro
+ *
+ * This class provides the base class for all cell based variables defined
+ * on the AMR mesh, accessible through the tree data structure.
+ * @sa ARTIO::HYDRO::data, ARTIO::HYDRO::empty_data
+ */
+template< typename TreeType_, typename ValueType_=double >
+class proto_data{
+
+public:
+	std::vector< std::vector<ValueType_> > m_var_array;
+protected:
+	TreeType_& m_tree;  //!< reference to underlying AMR tree structure
+	unsigned			m_cpu;			//!< the computational domain
+	unsigned
+		m_minlevel, 				//!< maximum refinement level to be read from file
+		m_maxlevel;					//!< minimum refinement level to be read from file
+	unsigned 			m_twotondim;//!< 2**ndim
+	unsigned			m_ilevel; 	//!< the refinement level
+
+	//! array holding the actual data
+
+public:
+	
+	//! constructor for the base class of all hydro data objects
+	/*! 
+	 * @param AMRtree reference to the underlying AMR tree data structure object
+	 */
+	explicit proto_data( TreeType_& AMRtree )
+	: m_tree(AMRtree), m_cpu( AMRtree.m_cpu ), 
+	  m_minlevel( AMRtree.m_minlevel ), m_maxlevel( AMRtree.m_maxlevel ),
+		m_twotondim( (unsigned)(pow(2, AMRtree.m_header.ndim)+0.5) )
+	{ }
+
+	//! access the value of the cells associated with the oct designated by the iterator
+	/*!
+	 * @param it the grid iterator pointing to the current oct
+	 * @param ind index of the child cell of the current oct (0..7)
+	 */
+	inline ValueType_& cell_value( const typename TreeType_::iterator& it, int ind )
+	{
+		unsigned ipos   = it.get_absolute_position();
+		unsigned ilevel = it.get_level();//-m_minlevel;
+		return (m_var_array[ilevel])[m_twotondim*ipos+ind];
+	}
+
+	//! access the value of the cells associated with the oct designated by the iterator
+	/*!
+	 * @param it the grid iterator pointing to the current oct
+	 * @param ind index of the child cell of the current oct (0..7)
+	 */
+	inline ValueType_& operator()( const typename TreeType_::iterator& it, int ind )
+	{	return cell_value(it,ind); }
+	
+	
+	//! combines all elements of this instance with that of another using a binary operator
+	/*!
+	 * @param o  the other data object with which to combine the elements
+	 * @param op the binary operator to be used in combining elements
+	 */
+	template<typename BinaryOperator_>
+	void combine_with( const proto_data<TreeType_,ValueType_>& o, const BinaryOperator_& op )
+	{
+		if( m_minlevel != o.m_minlevel || m_maxlevel != o.m_maxlevel 
+				|| m_twotondim != o.m_twotondim || &m_tree != &o.m_tree 
+				|| m_var_array.size() != o.m_var_array.size() ){
+			std::cerr << "this #levels=" << m_var_array.size() << ", other #levels=" << o.m_var_array.size() << std::endl;
+			throw std::runtime_error("Error: trying to combine incompatible mesh data.");
+		}
+		
+		for( unsigned ilvl=0; ilvl<m_var_array.size(); ++ilvl )
+		{
+			if( m_var_array[ilvl].size() != o.m_var_array[ilvl].size() ){
+				std::cerr << "ilvl=" << ilvl << ", this size=" << m_var_array[ilvl].size() << ", other size=" << o.m_var_array[ilvl].size() << std::endl;
+				throw std::runtime_error("Error: trying to combine incompatible mesh data.");
+			}
+			
+			for( unsigned i=0; i<m_var_array[ilvl].size(); ++i )
+				m_var_array[ilvl][i] = op( m_var_array[ilvl][i], o.m_var_array[ilvl][i] );
+		}
+		
+	}
+};
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+/*!
+ * @class ARTIO::HYDRO::empty_data
+ * @brief encapsulates additional (non-ARTIO) hydro data
+ *
+ * This class provides a wrapper for additional hydro variables created
+ * during post-processing. They are stored compatibly with the normal
+ * ARTIO hydro variables and thus can be accessed via the tree.
+ * @sa ARTIO::HYDRO::data, ARTIO::HYDRO::proto_data
+ */
+template< typename TreeType_, typename ValueType_=double >
+class empty_data : public proto_data<TreeType_,ValueType_>{
+	
+public:
+
+	//! constructor for an empty hydro variable defined on the AMR mesh
+	/*!
+	 * @param AMRtree the underlying hierarchical tree data structure
+	 * @param v the value with which the variable should be initialized, default is a (double) zero.
+	 */
+	explicit empty_data( TreeType_& AMRtree, ValueType_ v=(ValueType_)0.0 )
+	: proto_data<TreeType_,ValueType_>(AMRtree)
+	{
+		this->m_var_array.assign( this->m_maxlevel+1, std::vector<ValueType_>() );
+		
+		for( unsigned ilvl = this->m_minlevel; ilvl<=this->m_maxlevel; ++ilvl ){
+			typename TreeType_::iterator grid_it = this->m_tree.begin( ilvl );
+			while( grid_it != this->m_tree.end(ilvl) ){
+
+				for( unsigned j=0;j<this->m_twotondim;++j )
+					this->m_var_array[ilvl].push_back( v );
+
+				++grid_it;
+			}
+		}
+	}
+
+
+	//! write the new hydro variable to a ARTIO compatible output file
+	/*!
+	 * writes the hydro variable to a ARTIO compatible output file named
+	 * after the convention
+	 *  (path)/(basename)_(DOMAIN).out(snap_num)
+	 *
+	 * @param path the path where to store the files
+	 * @param basename the filename base string to prepend to the domain number
+	 * @param snap_num the number of the snapshot (default is zero).
+	 */
+	void save( std::string path, std::string basename, unsigned snap_num=0 )
+	{
+		char fullname[256];
+		sprintf(fullname,"%s/%s_%05d.out%05d",path.c_str(),basename.c_str(), snap_num, this->m_tree.m_cpu );
+		std::ofstream ofs( fullname, std::ios::binary|std::ios::trunc );
+		
+		
+		typename TreeType_::iterator it;
+		unsigned ncpu = this->m_tree.m_header.ncpu;
+		
+		//std::cerr << "ncpu = " << ncpu << std::endl;
+
+		for( unsigned ilvl = 0; ilvl<=this->m_maxlevel; ++ilvl ){
+			std::vector< std::vector<ValueType_> > 
+				temp1 (ncpu, std::vector<ValueType_>() ),
+				temp2 (ncpu, std::vector<ValueType_>() ),
+				temp3 (ncpu, std::vector<ValueType_>() ),
+				temp4 (ncpu, std::vector<ValueType_>() ),
+				temp5 (ncpu, std::vector<ValueType_>() ),
+				temp6 (ncpu, std::vector<ValueType_>() ),
+				temp7 (ncpu, std::vector<ValueType_>() ),
+				temp8 (ncpu, std::vector<ValueType_>() );
+			
+			it = this->m_tree.begin(ilvl);
+			
+			while( it!= this->m_tree.end(ilvl) )
+			{
+				temp1[ it.get_domain()-1 ].push_back( (*this)(it,0) );
+				temp2[ it.get_domain()-1 ].push_back( (*this)(it,1) );
+				temp3[ it.get_domain()-1 ].push_back( (*this)(it,2) );
+				temp4[ it.get_domain()-1 ].push_back( (*this)(it,3) );
+				temp5[ it.get_domain()-1 ].push_back( (*this)(it,4) );
+				temp6[ it.get_domain()-1 ].push_back( (*this)(it,5) );
+				temp7[ it.get_domain()-1 ].push_back( (*this)(it,6) );
+				temp8[ it.get_domain()-1 ].push_back( (*this)(it,7) );
+				++it;
+			}
+			
+			
+			for( unsigned icpu = 0; icpu<ncpu; ++icpu ){
+				unsigned nn;
+				
+				nn = temp1[icpu].size() * sizeof( ValueType_ );
+				
+				if( nn > 0 )
+				{
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp1[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp2[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp3[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp4[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp5[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp6[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp7[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+					ofs.write( (char*)&temp8[icpu][0], nn );
+					ofs.write( (char*)&nn, sizeof(unsigned) );
+				}
+			}
+		}
+	}
+	
+	//! reads an additional from a ARTIO compatible (single var) output file
+	/*!
+	 * @param basename the base string used for the variable
+	 */
+	void read( std::string basename )
+	{
+		char fullname[256];
+		sprintf(fullname,"%s_%05d.out%05d",basename.c_str(), this->m_tree.m_header.nout[0], this->m_tree.m_cpu );
+		std::ifstream ifs( fullname, std::ios::binary );
+
+		for( unsigned ilvl = this->m_minlevel; ilvl<=this->m_maxlevel; ++ilvl ){
+			unsigned nn = this->m_var_array[ilvl].size() * sizeof(ValueType_);
+			unsigned nnfile;
+			ifs.read( (char*)&nnfile, sizeof(unsigned) );
+			if( nn != nnfile ){
+				std::cerr << "Error: dimension mismatch between AMR tree and file data!" << std::endl;
+				std::cerr << "       found " << nnfile << ", expected " << nn << " in file \'" << fullname << "\'" << std::endl;
+				return;
+			}
+			ifs.read( (char*)&this->m_var_array[ilvl][0], nn );
+			ifs.read( (char*)&nn, sizeof(unsigned) );
+		}
+		ifs.close();
+	}
+
+};
+
+
+/**************************************************************************************\
+\**************************************************************************************/
+
+/*!
+ * @class data
+ * @brief encapsulates hydro data from a ARTIO simulation snapshot
+ *
+ * This class provides low-level read access to ARTIO hydro_XXXXX.out files.
+ * Data from a given list of computational domains can be read and is
+ * stored in internal datastructures.
+ * Access to cell position and threaded tree structure of the cell is provided
+ * through the member functions of class ARTIO_amr_level.
+ * @sa ARTIO_amr_level
+ */
+template< typename TreeType_, typename Real_=double >
+class data : public proto_data<TreeType_,Real_>{
+
+public:
+	struct header{
+		unsigned ncpu;		//!< number of CPUs in simulation
+		unsigned nvar;		//!< number of hydrodynamic variables
+		unsigned ndim;		//!< number of spatial dimensions
+		unsigned nlevelmax;	//!< maximum allowed refinement level
+		unsigned nboundary;	//!< number of boundary regions
+		double gamma;	//!< adiabatic exponent
+	};
+
+
+	std::string 	m_fname;		//!< the file name	
+	struct header	m_header;	 	//!< header meta data
+
+	const unsigned m_nvars; //!< number of variables stored in file
+
+	std::vector<std::string> m_varnames;	//!< names of the variables stored in file
+	std::map<std::string,unsigned> m_var_name_map; //!< a hash table for variable name to internal variable index
+
+
+protected:
+
+	//! generates a hydro_XXXX filename for specified cpu
+	std::string gen_fname( int icpu );
+
+	//! generate hydro_XXXXX filename from info filename
+	std::string rename_info2hydro( const std::string& info );
+
+	//! generate hydro_XXXXX filename from amr filename
+	std::string rename_amr2hydro( const std::string& info );
+
+	//! read header data containing meta information
+	void read_header( void );
+
+	//! get internal index for given variable string identifier
+	/*!
+	 * @param varname the string identifier of the hydro variable
+	 * @return internal variable index
+	 */
+	int get_var_idx( const std::string& varname )
+	{
+		int ivar;
+		std::map<std::string,unsigned>::iterator mit;
+		if( (mit=m_var_name_map.find(varname)) != m_var_name_map.end() )
+			ivar = (*mit).second;
+		else
+			throw std::runtime_error("ARTIO::HYDRO::data::get_var_idx :"\
+				" Error, cannot find variable named \'"+varname+"\'");
+		return ivar;
+	}
+
+	//! perform read operation of one hydro variable (internal use)
+	/*!
+	 * users should always call the read( std::string ) member function
+	 * and read variables through their string identifiers
+	 * @param var the index of the hydro variable
+	 */
+	void read( unsigned var );
+
+public:
+
+	//! constructor for hydro data
+	/*!
+	 * @param AMRtree reference to the underlying tree object
+	 */
+	explicit data( TreeType_& AMRtree )
+	: proto_data<TreeType_,Real_>( AMRtree ),
+	  m_fname( rename_amr2hydro(AMRtree.m_fname) ),
+		m_nvars( 6 )
+	{
+		read_header();
+
+		if( this->m_cpu > m_header.ncpu || this->m_cpu < 1 )
+			throw std::runtime_error("ARTIO::HYDRO::data : expect to read from out of range CPU.");
+
+		if( this->m_minlevel < 0 || this->m_maxlevel >= m_header.nlevelmax )
+			throw std::runtime_error("ARTIO::HYDRO::data : requested level is invalid.");
+
+		//m_twotondim = (unsigned)(pow(2,m_header.ndim)+0.5);
+
+		for( unsigned i=0; i<m_nvars; ++i ){
+                std::string tmp(artio_hydro_variables[i]);
+				m_var_name_map.insert( std::pair<std::string,unsigned>( tmp, i+1 ) );
+				m_varnames.push_back( tmp );
+			}
+	}
+
+	//! retrieve the names of variables available in the hydro data source
+	/*! Variables have an internal index but are accessed through their name, given as a string
+	 *  hydro data files may contain different variables depending on the kind of simulation run.
+	 * @param names An output iterator to which std::string designating available variables are sent
+	 * @return Final position of the output iterator
+	 */
+	template< typename _OutputIterator >