yt.geometry_tests / test_art_octre.py

from yt.mods import *
from art_header import read_header
from art_header import _count_art_octs
from art_header import _read_art_level_info

parameter_filename = 'ART/10MpcBox_csf512_a0.150.d'
root_grid_mask_offset,root_grid_offset,child_grid_offset,min_level,max_level = read_header(parameter_filename)
f = open(parameter_filename,'rb')
nhydrovars, inoll, level_oct_offsets, level_child_offsets =_count_art_octs(f, child_grid_offset,min_level, max_level)

NX = na.ones(3)*256
LE = na.array([0.0, 0.0, 0.0], dtype='float64')
RE = na.array([1.0, 1.0, 1.0], dtype='float64')

nocts = 256**3
pos = f.tell()
for level in range(1, max_level+1):
    # read level info into POS
    # 1 is the CPU, -1 is not used, 1 is the domain ID
    le, fln, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
    nocts += le.shape[0]
    print level, nLevel
print "NUMBER OF OCTS", nocts

art_oct = RAMSESOctreeContainer(NX, LE, RE)
art_oct.allocate_domains([nocts])

root_dx = (RE - LE) / NX
LL = LE + root_dx/2.0
RL = RE - root_dx/2.0
rpos = na.mgrid[ LL[0]:RL[0]:NX[0]*1j,
                 LL[1]:RL[1]:NX[1]*1j,
                 LL[2]:RL[2]:NX[2]*1j ]
rpos = na.vstack([p.ravel() for p in rpos]).T
print rpos.shape
art_oct.add(1, 0, -1, rpos, 1)

f.seek(pos)
for level in range(1, max_level):
    # read level info into POSx
    # 1 is the CPU, -1 is not used, 1 is the domain ID
    le, fln, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
    # int curdom, int curlevel, int ng,
    #       np.ndarray[np.float64_t, ndim=2] pos,
    #        int local_domain):
    fle = (le.astype("float64") / (NX*2**level)) * (RE - LE) + LE
    print "LEVEL", level, fle.min(), fle.max()
    art_oct.add(1, level, -1, fle, 1)

# Now let's test a selector
class Empty(object):
    pass

sp_obj = Empty()
sp_obj.center = [0.5, 0.5, 0.5]
sp_obj.radius = 0.1
sp = sel.SphereSelector(sp_obj)

mask = na.zeros(sp.nocts, dtype="bool")

cells = art_oct.count_cells(sp, mask)
nc = cells.sum()
print nc, sp.nocts
# 1 is domain_id
icoords = art_oct.icoords(1, mask, nc)
ires = art_oct.ires(1, mask, nc)

print icoords
print ires
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.