[Yt-svn] yt: Adding a reader to get the ART hierarchy relatively quickly....

[hg at spacepope.org]

hg Repository: yt
details:   yt/rev/80e00573f7a9
changeset: 3440:80e00573f7a9
user:      Matthew Turk <matthewturk at gmail.com>
date:
Wed Oct 13 16:38:35 2010 -0700
description:
Adding a reader to get the ART hierarchy relatively quickly.  This assumes some
things about endianness of the system you're on.

diffstat:

 yt/utilities/_amr_utils/endian_swap.h      |  14 +++++++
 yt/utilities/_amr_utils/fortran_reader.pyx |  61 ++++++++++++++++++++++++++++++
 2 files changed, 75 insertions(+), 0 deletions(-)

diffs (93 lines):

diff -r e061e70ed4ea -r 80e00573f7a9 yt/utilities/_amr_utils/endian_swap.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/endian_swap.h	Wed Oct 13 16:38:35 2010 -0700
@@ -0,0 +1,14 @@
+/* 
+   These macros are taken from Paul Bourke's page at
+   http://local.wasp.uwa.edu.au/~pbourke/dataformats/fortran/
+
+   The current year is 2010, and evidently we still have to deal with
+   endianness.
+*/
+
+#define SWAP_2(x) ( (((x) & 0xff) << 8) | ((unsigned short)(x) >> 8) )
+#define SWAP_4(x) ( ((x) << 24) | (((x) << 8) & 0x00ff0000) | \
+         (((x) >> 8) & 0x0000ff00) | ((x) >> 24) )
+#define FIX_SHORT(x) (*(unsigned short *)&(x) = SWAP_2(*(unsigned short *)&(x)))
+#define FIX_LONG(x) (*(unsigned *)&(x) = SWAP_4(*(unsigned *)&(x)))
+#define FIX_FLOAT(x) FIX_LONG(x)
diff -r e061e70ed4ea -r 80e00573f7a9 yt/utilities/_amr_utils/fortran_reader.pyx
--- a/yt/utilities/_amr_utils/fortran_reader.pyx	Wed Oct 13 16:02:27 2010 -0700
+++ b/yt/utilities/_amr_utils/fortran_reader.pyx	Wed Oct 13 16:38:35 2010 -0700
@@ -29,6 +29,14 @@
 
 from stdio cimport fopen, fclose, FILE
 
+cdef extern from "endian_swap.h":
+    void FIX_SHORT( unsigned short )
+    void FIX_LONG( unsigned )
+    void FIX_FLOAT( float )
+
+cdef extern from "alloca.h":
+    void *alloca(int)
+
 cdef extern from "stdio.h":
     cdef int SEEK_SET
     cdef int SEEK_CUR
@@ -68,3 +76,56 @@
             fread(<void *> (data + pos), 4, slab_size[0], f)
             pos += slab_size[0]
     return buffer
+
+def read_art_tree(char *fn, long offset,
+                  int min_level, int max_level,
+                  np.ndarray[np.int64_t, ndim=2] oct_indices,
+                  np.ndarray[np.int64_t, ndim=1] oct_levels,
+                  np.ndarray[np.int64_t, ndim=1] oct_parents
+                  ):
+    # This accepts the filename of the ART header and an integer offset that
+    # points to the start of the record *following* the reading of iOctFree and
+    # nOct.  For those following along at home, we only need to read:
+    #   iOctPr, iOctLv
+    cdef int nchild = 8
+    cdef int i, Lev, cell_ind, iOct, nLevel, nLevCells, ic1, next_record
+    cdef int iOctPs[3]
+    cdef int dummy_records[9]
+    cdef int readin
+    cdef FILE *f = fopen(fn, "rb")
+    fseek(f, offset, SEEK_SET)
+    cdef int * Level = <int *> alloca(sizeof(int)*(max_level-min_level+1))
+    cdef int * iNOLL = <int *> alloca(sizeof(int)*(max_level-min_level+1))
+    cdef int * iHOLL = <int *> alloca(sizeof(int)*(max_level-min_level+1))
+    cell_ind = 0
+    for Lev in xrange(min_level + 1, max_level + 1):
+        fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+        print "Reading Hierarchy for Level"
+        fread(&Level[Lev], sizeof(int), 1, f); FIX_LONG(Level[Lev])
+        fread(&iNOLL[Lev], sizeof(int), 1, f); FIX_LONG(iNOLL[Lev])
+        fread(&iHOLL[Lev], sizeof(int), 1, f); FIX_LONG(iHOLL[Lev])
+        fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+        iOct = iHOLL[Lev] - 1
+        nLevel = iNOLL[Lev]
+        for ic1 in xrange(nLevel):
+            #print readin, iOct, nLevel, sizeof(int) 
+            next_record = ftell(f)
+            fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+            next_record += readin + 2*sizeof(int)
+            fread(iOctPs, sizeof(int), 3, f)
+            FIX_LONG(iOctPs[0]); FIX_LONG(iOctPs[1]); FIX_LONG(iOctPs[2])
+            oct_indices[iOct, 0] = iOctPs[0]
+            oct_indices[iOct, 1] = iOctPs[1]
+            oct_indices[iOct, 2] = iOctPs[2]
+            fread(dummy_records, sizeof(int), 6, f) # skip Nb
+            fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+            oct_parents[iOct] = readin
+            fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+            oct_levels[iOct] = readin
+            fread(&iOct, sizeof(int), 1, f); FIX_LONG(iOct);
+            iOct -= 1
+            fseek(f, next_record, SEEK_SET)
+        fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+        next_record = (2*sizeof(int) + readin) * (nLevel * nchild)
+        next_record -= sizeof(int)
+        fseek(f, next_record, SEEK_CUR)