[Yt-svn] commit/yt: 2 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Sep 22 10:39:31 PDT 2011
2 new changesets in yt:
http://bitbucket.org/yt_analysis/yt/changeset/e78b2426cb13/
changeset: e78b2426cb13
branch: yt
user: MatthewTurk
date: 2011-09-22 19:38:57
summary: Adding FigureCanvasPS to _mpl_imports
affected #: 1 file (-1 bytes)
--- a/yt/visualization/_mpl_imports.py Sat Sep 03 11:48:19 2011 -0600
+++ b/yt/visualization/_mpl_imports.py Thu Sep 22 10:38:57 2011 -0700
@@ -17,6 +17,9 @@
from matplotlib.backends.backend_pdf import \
FigureCanvasPdf
+from matplotlib.backends.backend_ps import \
+ FigureCanvasPS
+
# Now we provide some convenience functions to get information about plots.
# With Matplotlib 0.98.x, the 'transforms' branch broke backwards
# compatibility. Despite that, the various packagers are plowing ahead with
http://bitbucket.org/yt_analysis/yt/changeset/41ed282c6e02/
changeset: 41ed282c6e02
branch: yt
user: MatthewTurk
date: 2011-09-22 19:39:25
summary: Merge
affected #: 18 files (-1 bytes)
--- a/CREDITS Thu Sep 22 10:38:57 2011 -0700
+++ b/CREDITS Thu Sep 22 10:39:25 2011 -0700
@@ -20,6 +20,7 @@
Casey Stark (caseywstark at gmail.com)
JC Passy (jcpassy at gmail.com)
Eve Lee (elee at cita.utoronto.ca)
+ Elizabeth Tasker (tasker at astro1.sci.hokudai.ac.jp)
We also include the Delaunay Triangulation module written by Robert Kern of
Enthought, the cmdln.py module by Trent Mick, and the progressbar module by
--- a/scripts/iyt Thu Sep 22 10:38:57 2011 -0700
+++ b/scripts/iyt Thu Sep 22 10:39:25 2011 -0700
@@ -3,6 +3,7 @@
from yt.mods import *
from yt.data_objects.data_containers import AMRData
namespace = locals().copy()
+namespace.pop("__builtins__", None)
doc = """\
--- a/yt/analysis_modules/halo_mass_function/halo_mass_function.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/analysis_modules/halo_mass_function/halo_mass_function.py Thu Sep 22 10:39:25 2011 -0700
@@ -700,7 +700,7 @@
def SQR(a):
return a*a
-def integrate_inf(fcn, error=1e-7, initial_guess=10):
+def integrate_inf(fcn, error=1e-3, initial_guess=10):
"""
Integrate a function *fcn* from zero to infinity, stopping when the answer
changes by less than *error*. Hopefully someday we can do something
--- a/yt/data_objects/data_containers.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/data_objects/data_containers.py Thu Sep 22 10:39:25 2011 -0700
@@ -40,7 +40,7 @@
from yt.data_objects.particle_io import particle_handler_registry
from yt.utilities.amr_utils import find_grids_in_inclined_box, \
grid_points_in_volume, planar_points_in_volume, VoxelTraversal, \
- QuadTree, get_box_grids_below_level
+ QuadTree, get_box_grids_below_level, ghost_zone_interpolate
from yt.utilities.data_point_utilities import CombineGrids, \
DataCubeRefine, DataCubeReplace, FillRegion, FillBuffer
from yt.utilities.definitions import axis_names, x_dict, y_dict
@@ -3104,6 +3104,9 @@
_type_name = "smoothed_covering_grid"
@wraps(AMRCoveringGridBase.__init__)
def __init__(self, *args, **kwargs):
+ self._base_dx = (
+ (self.pf.domain_right_edge - self.pf.domain_left_edge) /
+ self.pf.domain_dimensions.astype("float64"))
AMRCoveringGridBase.__init__(self, *args, **kwargs)
self._final_start_index = self.global_startindex
@@ -3155,8 +3158,10 @@
if self._use_pbar: pbar.finish()
def _update_level_state(self, level, field = None):
- dx = self.pf.h.select_grids(level)[0].dds
- for ax, v in zip('xyz', dx): self['cd%s'%ax] = v
+ dx = self._base_dx / self.pf.refine_by**level
+ self.data['cdx'] = dx[0]
+ self.data['cdy'] = dx[1]
+ self.data['cdz'] = dx[2]
LL = self.left_edge - self.pf.domain_left_edge
self._old_global_startindex = self.global_startindex
self.global_startindex = na.rint(LL / dx).astype('int64') - 1
@@ -3165,41 +3170,29 @@
if level == 0 and self.level > 0:
# We use one grid cell at LEAST, plus one buffer on all sides
idims = na.rint((self.right_edge-self.left_edge)/dx).astype('int64') + 2
- self[field] = na.zeros(idims,dtype='float64')-999
+ self.data[field] = na.zeros(idims,dtype='float64')-999
self._cur_dims = idims.astype("int32")
elif level == 0 and self.level == 0:
DLE = self.pf.domain_left_edge
self.global_startindex = na.array(na.floor(LL/ dx), dtype='int64')
idims = na.rint((self.right_edge-self.left_edge)/dx).astype('int64')
- self[field] = na.zeros(idims,dtype='float64')-999
+ self.data[field] = na.zeros(idims,dtype='float64')-999
self._cur_dims = idims.astype("int32")
def _refine(self, dlevel, field):
rf = float(self.pf.refine_by**dlevel)
- old_dims = na.array(self[field].shape) - 1
- old_left = (self._old_global_startindex + 0.5) * rf
- old_right = rf*old_dims + old_left
- old_bounds = [old_left[0], old_right[0],
- old_left[1], old_right[1],
- old_left[2], old_right[2]]
+ input_left = (self._old_global_startindex + 0.5) * rf
+ dx = na.fromiter((self['cd%s' % ax] for ax in 'xyz'), count=3, dtype='float64')
+ output_dims = na.rint((self.right_edge-self.left_edge)/dx).astype('int32') + 2
- dx = na.array([self['cd%s' % ax] for ax in 'xyz'], dtype='float64')
- new_dims = na.rint((self.right_edge-self.left_edge)/dx).astype('int64') + 2
+ self._cur_dims = output_dims
- # x, y, z are the new bounds
- x,y,z = (na.mgrid[0:new_dims[0], 0:new_dims[1], 0:new_dims[2]]
- ).astype('float64') + 0.5
- x += self.global_startindex[0]
- y += self.global_startindex[1]
- z += self.global_startindex[2]
- fake_grid = {'x':x,'y':y,'z':z}
-
- interpolator = TrilinearFieldInterpolator(
- self[field], old_bounds, ['x','y','z'],
- truncate = True)
- self._cur_dims = new_dims.astype("int32")
- self[field] = interpolator(fake_grid)
+ output_field = na.zeros(output_dims, dtype="float64")
+ output_left = self.global_startindex + 0.5
+ ghost_zone_interpolate(rf, self[field], input_left,
+ output_field, output_left)
+ self[field] = output_field
def _get_data_from_grid(self, grid, fields, level):
fields = ensure_list(fields)
--- a/yt/data_objects/grid_patch.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/data_objects/grid_patch.py Thu Sep 22 10:39:25 2011 -0700
@@ -71,7 +71,7 @@
level.
"""
- if self.start_index != None:
+ if self.start_index is not None:
return self.start_index
if self.Parent == None:
iLE = self.LeftEdge - self.pf.domain_left_edge
--- a/yt/data_objects/universal_fields.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/data_objects/universal_fields.py Thu Sep 22 10:39:25 2011 -0700
@@ -920,10 +920,9 @@
display_name=r"\mathrm{Particle}\/\mathrm{Density})")
def _MagneticEnergy(field,data):
- """WARNING WARNING WARNING: Units are not yet known to be
- correct. Trust the magnitude of this quantity at your own
- risk. However, it should just be a multiplicative offset from
- reality...
+ """This assumes that your front end has provided Bx, By, Bz in
+ units of Gauss. If you use MKS, make sure to write your own
+ MagneticEnergy field to deal with non-unitary \mu_0.
"""
return (data["Bx"]**2 + data["By"]**2 + data["Bz"]**2)/2.
add_field("MagneticEnergy",function=_MagneticEnergy,
--- a/yt/frontends/enzo/data_structures.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/frontends/enzo/data_structures.py Thu Sep 22 10:39:25 2011 -0700
@@ -220,10 +220,9 @@
# Sets are sorted, so that won't work!
def _parse_hierarchy(self):
def _next_token_line(token, f):
- line = f.readline()
- while token not in line:
- line = f.readline()
- return line.split()[2:]
+ for line in f:
+ if line.startswith(token):
+ return line.split()[2:]
if os.path.exists(self.hierarchy_filename[:-9] + "harrays"):
if self._parse_binary_hierarchy(): return
t1 = time.time()
@@ -234,7 +233,6 @@
self.grids[0].Level = 0
si, ei, LE, RE, fn, np = [], [], [], [], [], []
all = [si, ei, LE, RE, fn]
- f.readline() # Blank at top
pbar = get_pbar("Parsing Hierarchy", self.num_grids)
for grid_id in xrange(self.num_grids):
pbar.update(grid_id)
@@ -248,15 +246,11 @@
if nb > 0: fn[-1] = _next_token_line("BaryonFileName", f)
np.append(int(_next_token_line("NumberOfParticles", f)[0]))
if nb == 0 and np[-1] > 0: fn[-1] = _next_token_line("FileName", f)
- line = f.readline()
- while len(line) > 2:
+ for line in f:
+ if len(line) < 2: break
if line.startswith("Pointer:"):
vv = patt.findall(line)[0]
self.__pointer_handler(vv)
- line = f.readline()
- continue
- params = line.split()
- line = f.readline()
pbar.finish()
self._fill_arrays(ei, si, LE, RE, np)
self.grids = na.array(self.grids, dtype='object')
@@ -764,6 +758,8 @@
else:
if any("." in v or "e+" in v or "e-" in v for v in vals):
pcast = float
+ elif v == "inf":
+ pcast = str
else:
pcast = int
# Now we figure out what to do with it.
@@ -795,13 +791,18 @@
self.conversion_factors[data_labels[k]] = v
self.refine_by = self.parameters["RefineBy"]
self.dimensionality = self.parameters["TopGridRank"]
- self.domain_dimensions = self.parameters["TopGridDimensions"]
if self.dimensionality > 1:
- self.domain_left_edge = na.array(self.parameters["DomainLeftEdge"]).copy()
- self.domain_right_edge = na.array(self.parameters["DomainRightEdge"]).copy()
+ self.domain_dimensions = self.parameters["TopGridDimensions"]
+ self.domain_left_edge = na.array(self.parameters["DomainLeftEdge"],
+ "float64").copy()
+ self.domain_right_edge = na.array(self.parameters["DomainRightEdge"],
+ "float64").copy()
else:
- self.domain_left_edge = na.array(self.parameters["DomainLeftEdge"])
- self.domain_right_edge = na.array(self.parameters["DomainRightEdge"])
+ self.domain_left_edge = na.array(self.parameters["DomainLeftEdge"],
+ "float64")
+ self.domain_right_edge = na.array(self.parameters["DomainRightEdge"],
+ "float64")
+ self.domain_dimensions = na.array([self.parameters["TopGridDimensions"],1,1])
self.current_time = self.parameters["InitialTime"]
# To be enabled when we can break old pickles:
--- a/yt/frontends/enzo/fields.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/frontends/enzo/fields.py Thu Sep 22 10:39:25 2011 -0700
@@ -380,7 +380,7 @@
for field in ['Bx','By','Bz']:
f = EnzoFieldInfo[field]
f._convert_function=_convertBfield
- f._units=r"\mathrm{Gau\ss}"
+ f._units=r"\mathrm{Gauss}"
f.take_log=False
def _Bmag(field, data):
@@ -388,7 +388,7 @@
"""
return na.sqrt(data['Bx']**2 + data['By']**2 + data['Bz']**2)
-add_field("Bmag", function=_Bmag,display_name=r"|B|",units=r"\mathrm{Gau\ss}")
+add_field("Bmag", function=_Bmag,display_name=r"|B|",units=r"\mathrm{Gauss}")
#
@@ -476,7 +476,7 @@
for field in ['Bx','By','Bz']:
f = EnzoFieldInfo[field]
f._convert_function=_convertBfield
- f._units=r"\mathrm{Gau\ss}"
+ f._units=r"\mathrm{Gauss}"
f.take_log=False
def _Bmag(field, data):
@@ -484,4 +484,4 @@
"""
return na.sqrt(data['Bx']**2 + data['By']**2 + data['Bz']**2)
-add_field("Bmag", function=_Bmag,display_name=r"|B|",units=r"\mathrm{Gau\ss}")
+add_field("Bmag", function=_Bmag,display_name=r"|B|",units=r"\mathrm{Gauss}")
--- a/yt/utilities/_amr_utils/FixedInterpolator.c Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/_amr_utils/FixedInterpolator.c Thu Sep 22 10:39:25 2011 -0700
@@ -85,36 +85,21 @@
void vertex_interp(npy_float64 v1, npy_float64 v2, npy_float64 isovalue,
npy_float64 vl[3], npy_float64 dds[3],
npy_float64 x, npy_float64 y, npy_float64 z,
- int x0, int y0, int z0, int direction)
+ int vind1, int vind2)
{
/*if (fabs(isovalue - v1) < 0.000001) return 0.0;
if (fabs(isovalue - v2) < 0.000001) return 1.0;
if (fabs(v1 - v2) < 0.000001) return 0.0;*/
- npy_float64 mu = (isovalue - v1) / (v2 - v1);
+ int i;
+ static npy_float64 cverts[8][3] =
+ {{0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
+ {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}};
+
+ npy_float64 mu = ((isovalue - v1) / (v2 - v1));
vl[0] = x; vl[1] = y; vl[2] = z;
- if (x0 == 1) vl[0] += dds[0];
- if (y0 == 1) vl[1] += dds[1];
- if (z0 == 1) vl[2] += dds[2];
- switch (direction) {
- case -1:
- vl[0] -= (1.0 - mu) * dds[0];
- break;
- case 1:
- vl[0] += mu * dds[0];
- break;
- case -2:
- vl[1] -= (1.0 - mu) * dds[1];
- break;
- case 2:
- vl[1] += mu * dds[1];
- break;
- case -3:
- vl[2] -= (1.0 - mu) * dds[2];
- break;
- case 3:
- vl[2] += mu * dds[2];
- break;
- }
+ for (i=0;i<3;i++)
+ vl[i] += dds[i] * cverts[vind1][i]
+ + dds[i] * mu*(cverts[vind2][i] - cverts[vind1][i]);
}
npy_float64 trilinear_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
--- a/yt/utilities/_amr_utils/FixedInterpolator.h Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/_amr_utils/FixedInterpolator.h Thu Sep 22 10:39:25 2011 -0700
@@ -47,4 +47,4 @@
void vertex_interp(npy_float64 v1, npy_float64 v2, npy_float64 isovalue,
npy_float64 vl[3], npy_float64 dds[3],
npy_float64 x, npy_float64 y, npy_float64 z,
- int x0, int y0, int z0, int direction);
+ int vind1, int vind2);
--- a/yt/utilities/_amr_utils/Interpolators.pyx Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/_amr_utils/Interpolators.pyx Thu Sep 22 10:39:25 2011 -0700
@@ -27,6 +27,8 @@
cimport numpy as np
cimport cython
+ at cython.cdivision(True)
+ at cython.wraparound(False)
@cython.boundscheck(False)
def UnilinearlyInterpolate(np.ndarray[np.float64_t, ndim=1] table,
np.ndarray[np.float64_t, ndim=1] x_vals,
@@ -44,6 +46,8 @@
output[i] = table[x_i ] * (xm) \
+ table[x_i+1] * (xp)
+ at cython.cdivision(True)
+ at cython.wraparound(False)
@cython.boundscheck(False)
def BilinearlyInterpolate(np.ndarray[np.float64_t, ndim=2] table,
np.ndarray[np.float64_t, ndim=1] x_vals,
@@ -73,6 +77,8 @@
+ table[x_i , y_i+1] * (xm*yp) \
+ table[x_i+1, y_i+1] * (xp*yp)
+ at cython.cdivision(True)
+ at cython.wraparound(False)
@cython.boundscheck(False)
def TrilinearlyInterpolate(np.ndarray[np.float64_t, ndim=3] table,
np.ndarray[np.float64_t, ndim=1] x_vals,
@@ -114,3 +120,54 @@
+ table[x_i ,y_i+1,z_i+1] * (xm*yp*zp) \
+ table[x_i+1,y_i+1,z_i ] * (xp*yp*zm) \
+ table[x_i+1,y_i+1,z_i+1] * (xp*yp*zp)
+
+ at cython.cdivision(True)
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def ghost_zone_interpolate(int rf,
+ np.ndarray[np.float64_t, ndim=3] input_field,
+ np.ndarray[np.float64_t, ndim=1] input_left,
+ np.ndarray[np.float64_t, ndim=3] output_field,
+ np.ndarray[np.float64_t, ndim=1] output_left):
+ cdef int oi, oj, ok
+ cdef int ii, ij, ik
+ cdef np.float64_t xp, xm, yp, ym, zp, zm
+ cdef np.float64_t ods[3], ids[3], iids[3]
+ cdef np.float64_t opos[3], ropos[3], temp
+ cdef int i, j
+ for i in range(3):
+ temp = input_left[i] + (rf * (input_field.shape[i] - 1))
+ ids[i] = (temp - input_left[i])/(input_field.shape[i]-1)
+ temp = output_left[i] + output_field.shape[i] - 1
+ ods[i] = (temp - output_left[i])/(output_field.shape[i]-1)
+ iids[i] = 1.0/ids[i]
+ opos[0] = output_left[0]
+ for oi in range(output_field.shape[0]):
+ ropos[0] = ((opos[0] - input_left[0]) * iids[0])
+ ii = iclip(<int> ropos[0], 0, input_field.shape[0] - 2)
+ xp = ropos[0] - ii
+ xm = 1.0 - xp
+ opos[1] = output_left[1]
+ for oj in range(output_field.shape[1]):
+ ropos[1] = ((opos[1] - input_left[1]) * iids[1])
+ ij = iclip(<int> ropos[1], 0, input_field.shape[1] - 2)
+ yp = ropos[1] - ij
+ ym = 1.0 - yp
+ opos[2] = output_left[2]
+ for ok in range(output_field.shape[2]):
+ ropos[2] = ((opos[2] - input_left[2]) * iids[2])
+ ik = iclip(<int> ropos[2], 0, input_field.shape[2] - 2)
+ zp = ropos[2] - ik
+ zm = 1.0 - zp
+ output_field[oi,oj,ok] = \
+ input_field[ii ,ij ,ik ] * (xm*ym*zm) \
+ + input_field[ii+1,ij ,ik ] * (xp*ym*zm) \
+ + input_field[ii ,ij+1,ik ] * (xm*yp*zm) \
+ + input_field[ii ,ij ,ik+1] * (xm*ym*zp) \
+ + input_field[ii+1,ij ,ik+1] * (xp*ym*zp) \
+ + input_field[ii ,ij+1,ik+1] * (xm*yp*zp) \
+ + input_field[ii+1,ij+1,ik ] * (xp*yp*zm) \
+ + input_field[ii+1,ij+1,ik+1] * (xp*yp*zp)
+ opos[2] += ods[2]
+ opos[1] += ods[1]
+ opos[0] += ods[0]
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Sep 22 10:39:25 2011 -0700
@@ -120,7 +120,7 @@
void vertex_interp(np.float64_t v1, np.float64_t v2, np.float64_t isovalue,
np.float64_t vl[3], np.float64_t dds[3],
np.float64_t x, np.float64_t y, np.float64_t z,
- int x0, int y0, int z0, int direction)
+ int vind1, int vind2)
#cdef extern int *edge_table
#cdef extern int **tri_table
@@ -1138,6 +1138,7 @@
[0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]]
+
cdef int i, j, k, n
cdef int offset
cdef np.float64_t gv[8]
@@ -1151,64 +1152,57 @@
cdef Triangle *current = NULL
x = self.left_edge[0]
for i in range(self.dims[0]):
- x += self.dds[0]
y = self.left_edge[1]
for j in range(self.dims[1]):
- y += self.dds[1]
z = self.left_edge[2]
for k in range(self.dims[2]):
- z += self.dds[2]
offset = i * (self.dims[1] + 1) * (self.dims[2] + 1) \
+ j * (self.dims[2] + 1) + k
intdata = self.data[field_id] + offset
offset_fill(self.dims, intdata, gv)
cubeindex = 0
- if gv[0] < isovalue: cubeindex |= 1
- if gv[1] < isovalue: cubeindex |= 2
- if gv[2] < isovalue: cubeindex |= 4
- if gv[3] < isovalue: cubeindex |= 8
- if gv[4] < isovalue: cubeindex |= 16
- if gv[5] < isovalue: cubeindex |= 32
- if gv[6] < isovalue: cubeindex |= 64
- if gv[7] < isovalue: cubeindex |= 128
- #print cubeindex
- if edge_table[cubeindex] == 0: continue
+ for n in range(8):
+ if gv[n] < isovalue:
+ cubeindex |= (1 << n)
+ if edge_table[cubeindex] == 0:
+ z += self.dds[2]
+ continue
if (edge_table[cubeindex] & 1): # 0,0,0 with 1,0,0
vertex_interp(gv[0], gv[1], isovalue, vertlist[0],
- self.dds, x, y, z, 0,0,0,1)
+ self.dds, x, y, z, 0, 1)
if (edge_table[cubeindex] & 2): # 1,0,0 with 1,1,0
vertex_interp(gv[1], gv[2], isovalue, vertlist[1],
- self.dds, x, y, z, 1,0,0,2)
+ self.dds, x, y, z, 1, 2)
if (edge_table[cubeindex] & 4): # 1,1,0 with 0,1,0
vertex_interp(gv[2], gv[3], isovalue, vertlist[2],
- self.dds, x, y, z, 1,0,0,-1)
+ self.dds, x, y, z, 2, 3)
if (edge_table[cubeindex] & 8): # 0,1,0 with 0,0,0
vertex_interp(gv[3], gv[0], isovalue, vertlist[3],
- self.dds, x, y, z, 0,1,0,-2)
+ self.dds, x, y, z, 3, 0)
if (edge_table[cubeindex] & 16): # 0,0,1 with 1,0,1
vertex_interp(gv[4], gv[5], isovalue, vertlist[4],
- self.dds, x, y, z, 0,0,1,1)
+ self.dds, x, y, z, 4, 5)
if (edge_table[cubeindex] & 32): # 1,0,1 with 1,1,1
vertex_interp(gv[5], gv[6], isovalue, vertlist[5],
- self.dds, x, y, z, 1,0,1,2)
+ self.dds, x, y, z, 5, 6)
if (edge_table[cubeindex] & 64): # 1,1,1 with 0,1,1
vertex_interp(gv[6], gv[7], isovalue, vertlist[6],
- self.dds, x, y, z, 1,1,1,-1)
+ self.dds, x, y, z, 6, 7)
if (edge_table[cubeindex] & 128): # 0,1,1 with 0,0,1
vertex_interp(gv[7], gv[4], isovalue, vertlist[7],
- self.dds, x, y, z, 0,1,1,-2)
+ self.dds, x, y, z, 7, 4)
if (edge_table[cubeindex] & 256): # 0,0,0 with 0,0,1
vertex_interp(gv[0], gv[4], isovalue, vertlist[8],
- self.dds, x, y, z, 0,0,0,3)
+ self.dds, x, y, z, 0, 4)
if (edge_table[cubeindex] & 512): # 1,0,0 with 1,0,1
vertex_interp(gv[1], gv[5], isovalue, vertlist[9],
- self.dds, x, y, z, 1,0,0,3)
+ self.dds, x, y, z, 1, 5)
if (edge_table[cubeindex] & 1024): # 1,1,0 with 1,1,1
vertex_interp(gv[2], gv[6], isovalue, vertlist[10],
- self.dds, x, y, z, 1,1,0,3)
+ self.dds, x, y, z, 2, 6)
if (edge_table[cubeindex] & 2048): # 0,1,0 with 0,1,1
vertex_interp(gv[3], gv[7], isovalue, vertlist[11],
- self.dds, x, y, z, 0,1,0,3)
+ self.dds, x, y, z, 3, 7)
n = 0
while 1:
current = AddTriangle(current,
@@ -1219,6 +1213,9 @@
if first == NULL: first = current
n += 3
if tri_table[cubeindex][n] == -1: break
+ z += self.dds[2]
+ y += self.dds[1]
+ x += self.dds[0]
# Hallo, we are all done.
cdef np.ndarray[np.float64_t, ndim=2] vertices
vertices = np.zeros((ntriang*3,3), dtype='float64')
--- a/yt/utilities/amr_kdtree/amr_kdtree.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/amr_kdtree/amr_kdtree.py Thu Sep 22 10:39:25 2011 -0700
@@ -60,6 +60,38 @@
def _rchild_id(id): return (id<<1) + 2
def _parent_id(id): return (id-1)>>1
+steps = na.array([[-1, -1, -1],
+ [-1, -1, 0],
+ [-1, -1, 1],
+ [-1, 0, -1],
+ [-1, 0, 0],
+ [-1, 0, 1],
+ [-1, 1, -1],
+ [-1, 1, 0],
+ [-1, 1, 1],
+
+ [ 0, -1, -1],
+ [ 0, -1, 0],
+ [ 0, -1, 1],
+ [ 0, 0, -1],
+ # [ 0, 0, 0],
+ [ 0, 0, 1],
+ [ 0, 1, -1],
+ [ 0, 1, 0],
+ [ 0, 1, 1],
+
+ [ 1, -1, -1],
+ [ 1, -1, 0],
+ [ 1, -1, 1],
+ [ 1, 0, -1],
+ [ 1, 0, 0],
+ [ 1, 0, 1],
+ [ 1, 1, -1],
+ [ 1, 1, 0],
+ [ 1, 1, 1]
+ ])
+
+
class MasterNode(object):
r"""
A MasterNode object is the building block of the AMR kd-Tree.
@@ -259,6 +291,7 @@
self.current_split_dim = 0
self.pf = pf
+ self.sdx = self.pf.h.get_smallest_dx()
self._id_offset = pf.h.grids[0]._id_offset
if nprocs > len(pf.h.grids):
mylog.info('Parallel rendering requires that the number of \n \
@@ -487,9 +520,97 @@
brick, le=le, re=re, periodic=periodic)]
return grids
+ def locate_neighbors_from_position(self, position):
+ r"""Given a position, finds the 26 neighbor grids
+ and cell indices.
+
+ This is a mostly a wrapper for locate_neighbors.
+
+ Parameters
+ ----------
+ position: array-like
+ Position of interest
+
+ Returns
+ -------
+ grids: Numpy array of Grid objects
+ cis: List of neighbor cell index tuples
+
+ Both of these are neighbors that, relative to the current cell
+ index (i,j,k), are ordered as:
+
+ (i-1, j-1, k-1), (i-1, j-1, k ), (i-1, j-1, k+1), ...
+ (i-1, j , k-1), (i-1, j , k ), (i-1, j , k+1), ...
+ (i+1, j+1, k-1), (i-1, j-1, k ), (i+1, j+1, k+1)
+
+ That is they start from the lower left and proceed to upper
+ right varying the third index most frequently. Note that the
+ center cell (i,j,k) is ommitted.
+
+ """
+ position = na.array(position)
+ grid = self.locate_brick(position).grid
+ ci = ((position-grid.LeftEdge)/grid.dds).astype('int64')
+ return self.locate_neighbors(grid,ci)
+
+ def locate_neighbors(self, grid, ci):
+ r"""Given a grid and cell index, finds the 26 neighbor grids
+ and cell indices.
+
+ Parameters
+ ----------
+ grid: Grid Object
+ Grid containing the cell of interest
+ ci: array-like
+ The cell index of the cell of interest
+
+ Returns
+ -------
+ grids: Numpy array of Grid objects
+ cis: List of neighbor cell index tuples
+
+ Both of these are neighbors that, relative to the current cell
+ index (i,j,k), are ordered as:
+
+ (i-1, j-1, k-1), (i-1, j-1, k ), (i-1, j-1, k+1), ...
+ (i-1, j , k-1), (i-1, j , k ), (i-1, j , k+1), ...
+ (i+1, j+1, k-1), (i-1, j-1, k ), (i+1, j+1, k+1)
+
+ That is they start from the lower left and proceed to upper
+ right varying the third index most frequently. Note that the
+ center cell (i,j,k) is ommitted.
+
+ """
+ ci = na.array(ci)
+ center_dds = grid.dds
+ position = grid.LeftEdge + (na.array(ci)+0.5)*grid.dds
+ grids = na.empty(26, dtype='object')
+ cis = na.empty([26,3], dtype='int64')
+ offs = 0.5*(center_dds + self.sdx)
+
+ new_cis = ci + steps
+ in_grid = na.all((new_cis >=0)*
+ (new_cis < grid.ActiveDimensions),axis=1)
+ new_positions = position + steps*offs
+ grids[in_grid] = grid
+
+ get_them = na.argwhere(in_grid != True).ravel()
+ cis[in_grid] = new_cis[in_grid]
+
+ if (in_grid != True).sum()>0:
+ grids[in_grid != True] = \
+ [self.locate_brick(new_positions[i]).grid for i in get_them]
+ cis[in_grid != True] = \
+ [(new_positions[i]-grids[i].LeftEdge)/
+ grids[i].dds for i in get_them]
+ cis = [tuple(ci) for ci in cis]
+ return grids, cis
+
def locate_brick(self, position):
r"""Given a position, find the node that contains it.
+ Will modify the position to account for periodicity.
+
Parameters
----------
pos: array_like
@@ -502,11 +623,17 @@
"""
node = self.tree
+ w = self.pf.domain_right_edge-self.pf.domain_left_edge
+ for i in range(3):
+ if position[i] < self.pf.domain_left_edge[i]:
+ position[i] += w[i]
+ if position[i] > self.pf.domain_right_edge[i]:
+ position[i] -= w[i]
while True:
if node.grid is not None:
return node
else:
- if position[node.split_ax] <= node.split_pos:
+ if position[node.split_ax] < node.split_pos:
node = node.left_child
else:
node = node.right_child
@@ -547,11 +674,14 @@
data = [d[current_node.li[0]:current_node.ri[0]+1,
current_node.li[1]:current_node.ri[1]+1,
current_node.li[2]:current_node.ri[2]+1].copy() for d in dds]
-
- current_node.brick = PartitionedGrid(current_node.grid.id, len(self.fields), data,
- current_node.l_corner.copy(),
- current_node.r_corner.copy(),
- current_node.dims.astype('int64'))
+
+ if na.any(current_node.r_corner-current_node.l_corner == 0):
+ current_node.brick = None
+ else:
+ current_node.brick = PartitionedGrid(current_node.grid.id, len(self.fields), data,
+ current_node.l_corner.copy(),
+ current_node.r_corner.copy(),
+ current_node.dims.astype('int64'))
self.bricks.append(current_node.brick)
self.brick_dimensions.append(current_node.dims)
self.bricks = na.array(self.bricks)
@@ -1010,7 +1140,8 @@
for node in self.viewpoint_traverse(viewpoint):
if node.grid is not None:
- yield node.brick
+ if node.brick is not None:
+ yield node.brick
self.reduce_tree_images(self.tree, front_center)
self._barrier()
--- a/yt/utilities/command_line.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/utilities/command_line.py Thu Sep 22 10:39:25 2011 -0700
@@ -173,9 +173,9 @@
action="store_true",
dest="enhance", default=False,
help="Enhance!"),
- range = dict(short="", long="--range",
+ valrange = dict(short="-r", long="--range",
action="store", type="float",
- dest="range", default=None,
+ dest="valrange", default=None,
nargs=2,
help="Range, space separated"),
up = dict(short="", long="--up",
@@ -1481,7 +1481,7 @@
@add_cmd_options(["width", "unit", "center","enhance",'outputfn',
"field", "cmap", "contours", "viewpoint",
- "pixels","up","range","log","contour_width"])
+ "pixels","up","valrange","log","contour_width"])
@check_args
def do_render(self, subcmd, opts, arg):
"""
@@ -1528,13 +1528,14 @@
if log is None:
log = True
- if opts.range is None:
+ myrange = opts.valrange
+ if myrange is None:
roi = pf.h.region(center, center-width, center+width)
mi, ma = roi.quantities['Extrema'](field)[0]
if log:
mi, ma = na.log10(mi), na.log10(ma)
else:
- mi, ma = range[0], range[1]
+ mi, ma = myrange[0], myrange[1]
n_contours = opts.contours
if n_contours is None:
@@ -1550,7 +1551,7 @@
cam = pf.h.camera(center, L, width, (N,N), transfer_function=tf)
image = cam.snapshot()
-
+
if opts.enhance:
for i in range(3):
image[:,:,i] = image[:,:,i]/(image[:,:,i].mean() + 5.*image[:,:,i].std())
@@ -1561,7 +1562,9 @@
save_name = "%s"%pf+"_"+field+"_rendering.png"
if not '.png' in save_name:
save_name += '.png'
- write_bitmap(image,save_name)
+ if cam._mpi_get_rank() != -1:
+ write_bitmap(image,save_name)
+
def run_main():
for co in ["--parallel", "--paste"]:
--- a/yt/visualization/color_maps.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/visualization/color_maps.py Thu Sep 22 10:39:25 2011 -0700
@@ -105,6 +105,25 @@
add_cmap('black_green', cdict)
+# This one comes from
+# http://permalink.gmane.org/gmane.comp.python.matplotlib.devel/10518
+# and is an implementation of http://arxiv.org/abs/1108.5083
+#
+
+# cubehelix parameters
+_gamma_cubehelix = 1.0
+_s_cubehelix = 0.5
+_r_cubehelix = -1.5
+_h_cubehelix = 1.0
+
+_cubehelix_data = {
+ 'red': lambda x: x**_gamma_cubehelix + (_h_cubehelix * x**_gamma_cubehelix * (1 - x**_gamma_cubehelix) / 2) * (-0.14861 * na.cos(2 * na.pi * (_s_cubehelix / 3 + _r_cubehelix * x)) + 1.78277 * na.sin(2 * na.pi * (_s_cubehelix / 3 + _r_cubehelix * x))),
+ 'green': lambda x: x**_gamma_cubehelix + (_h_cubehelix * x**_gamma_cubehelix * (1 - x**_gamma_cubehelix) / 2) * (-0.29227 * na.cos(2 * na.pi * (_s_cubehelix / 3 + _r_cubehelix * x)) - 0.90649 * na.sin(2 * na.pi * (_s_cubehelix / 3 + _r_cubehelix * x))),
+ 'blue': lambda x: x**_gamma_cubehelix + (_h_cubehelix * x**_gamma_cubehelix * (1 - x**_gamma_cubehelix) / 2) * (1.97294 * na.cos(2 * na.pi * (_s_cubehelix / 3 + _r_cubehelix * x))),
+}
+
+add_cmap("cubehelix", _cubehelix_data)
+
def _extract_lookup_table(cmap_name):
cmap = mcm.get_cmap(cmap_name)
if not cmap._isinit: cmap._init()
--- a/yt/visualization/eps_writer.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/visualization/eps_writer.py Thu Sep 22 10:39:25 2011 -0700
@@ -719,7 +719,7 @@
def multiplot(ncol, nrow, yt_plots=None, images=None, xranges=None,
yranges=None, xlabels=None, ylabels=None, colorbars=None,
shrink_cb=0.95, figsize=(8,8), margins=(0,0), titles=None,
- savefig=None, yt_nocbar=False, bare_axes=False,
+ savefig=None, format="eps", yt_nocbar=False, bare_axes=False,
cb_flags=None):
r"""Convenience routine to create a multi-panel figure from yt plots or
JPEGs. The images are first placed from the origin, and then
@@ -756,6 +756,8 @@
Titles that are placed in textboxes in each panel.
savefig : string
Name of the saved file without the extension.
+ format : string
+ File format of the figure. eps or pdf accepted.
yt_nocbar : boolean
Flag to indicate whether or not colorbars are created.
bare_axes : boolean
@@ -908,7 +910,7 @@
shrink=shrink_cb)
if savefig != None:
- d.save_fig(savefig)
+ d.save_fig(savefig, format=format)
return d
--- a/yt/visualization/plot_modifications.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/visualization/plot_modifications.py Thu Sep 22 10:39:25 2011 -0700
@@ -60,17 +60,23 @@
class VelocityCallback(PlotCallback):
_type_name = "velocity"
- def __init__(self, factor=16, scale=None, scale_units=None):
+ def __init__(self, factor=16, scale=None, scale_units=None, normalize=False):
"""
Adds a 'quiver' plot of velocity to the plot, skipping all but
- every *factor* datapoint
- *scale* is the data units per arrow length unit using *scale_units*
- (see matplotlib.axes.Axes.quiver for more info)
+ every *factor* datapoint. *scale* is the data units per arrow
+ length unit using *scale_units* (see
+ matplotlib.axes.Axes.quiver for more info). if *normalize* is
+ True, the velocity fields will be scaled by their local
+ (in-plane) length, allowing morphological features to be more
+ clearly seen for fields with substantial variation in field
+ strength (normalize is not implemented and thus ignored for
+ Cutting Planes).
"""
PlotCallback.__init__(self)
self.factor = factor
self.scale = scale
self.scale_units = scale_units
+ self.normalize = normalize
def __call__(self, plot):
# Instantiation of these is cheap
@@ -81,20 +87,26 @@
else:
xv = "%s-velocity" % (x_names[plot.data.axis])
yv = "%s-velocity" % (y_names[plot.data.axis])
- qcb = QuiverCallback(xv, yv, self.factor, scale=self.scale, scale_units=self.scale_units)
+ qcb = QuiverCallback(xv, yv, self.factor, scale=self.scale, scale_units=self.scale_units, normalize=self.normalize)
return qcb(plot)
class MagFieldCallback(PlotCallback):
_type_name = "magnetic_field"
- def __init__(self, factor=16, scale=None, scale_units=None):
+ def __init__(self, factor=16, scale=None, scale_units=None, normalize=False):
"""
Adds a 'quiver' plot of magnetic field to the plot, skipping all but
- every *factor* datapoint
+ every *factor* datapoint. *scale* is the data units per arrow
+ length unit using *scale_units* (see
+ matplotlib.axes.Axes.quiver for more info). if *normalize* is
+ True, the magnetic fields will be scaled by their local
+ (in-plane) length, allowing morphological features to be more
+ clearly seen for fields with substantial variation in field strength.
"""
PlotCallback.__init__(self)
self.factor = factor
self.scale = scale
self.scale_units = scale_units
+ self.normalize = normalize
def __call__(self, plot):
# Instantiation of these is cheap
@@ -103,12 +115,12 @@
else:
xv = "B%s" % (x_names[plot.data.axis])
yv = "B%s" % (y_names[plot.data.axis])
- qcb = QuiverCallback(xv, yv, self.factor, scale=self.scale, scale_units=self.scale_units)
+ qcb = QuiverCallback(xv, yv, self.factor, scale=self.scale, scale_units=self.scale_units, normalize=self.normalize)
return qcb(plot)
class QuiverCallback(PlotCallback):
_type_name = "quiver"
- def __init__(self, field_x, field_y, factor, scale=None, scale_units=None):
+ def __init__(self, field_x, field_y, factor, scale=None, scale_units=None, normalize=False):
"""
Adds a 'quiver' plot to any plot, using the *field_x* and *field_y*
from the associated data, skipping every *factor* datapoints
@@ -122,6 +134,7 @@
self.factor = factor
self.scale = scale
self.scale_units = scale_units
+ self.normalize = normalize
def __call__(self, plot):
x0, x1 = plot.xlim
@@ -147,6 +160,10 @@
(x0, x1, y0, y1),).transpose()
X = na.mgrid[0:plot.image._A.shape[0]-1:nx*1j]# + 0.5*factor
Y = na.mgrid[0:plot.image._A.shape[1]-1:ny*1j]# + 0.5*factor
+ if self.normalize:
+ nn = na.sqrt(pixX**2 + pixY**2)
+ pixX /= nn
+ pixY /= nn
plot._axes.quiver(X,Y, pixX, pixY, scale=self.scale, scale_units=self.scale_units)
plot._axes.set_xlim(xx0,xx1)
plot._axes.set_ylim(yy0,yy1)
--- a/yt/visualization/volume_rendering/grid_partitioner.py Thu Sep 22 10:38:57 2011 -0700
+++ b/yt/visualization/volume_rendering/grid_partitioner.py Thu Sep 22 10:39:25 2011 -0700
@@ -195,7 +195,7 @@
def initialize_source(self):
pass
- def traverse(self, back, front):
+ def traverse(self, back, front, image):
for b in self.bricks: yield b
def reset_cast(self):
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list