[Yt-svn] yt-commit r371 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Tue Jan 29 16:38:12 PST 2008
Author: mturk
Date: Tue Jan 29 16:38:11 2008
New Revision: 371
URL: http://yt.spacepope.org/changeset/371
Log:
Added a new cylinder datatype, which closes #55. Note that I am slightly
inconsistent -- it's a cylinder datatype, but it uses the .disk(...) method for
instantiation.
Fixed a print statement in the profiles.
Changed cut masks to be more conservative with memory by using just boolean
arrays rather than indices, and then un-did all of that good by adding an
automagical caching mechanism, which I will soon put in a parameter to shut
off.
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/DerivedFields.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/Profiles.py
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Tue Jan 29 16:38:11 2008
@@ -34,6 +34,14 @@
return tr
return save_state
+def cache_mask(func):
+ def check_cache(self, grid):
+ if not self._cut_masks.has_key(grid.id):
+ cm = func(self, grid)
+ self._cut_masks[grid.id] = cm
+ return self._cut_masks[grid.id]
+ return check_cache
+
class EnzoData:
"""
@@ -492,6 +500,7 @@
# Recall that the projection of the distance vector from a point
# onto the normal vector of a plane is:
# D = (a x_0 + b y_0 + c z_0 + d)/sqrt(a^2+b^2+c^2)
+ # @todo: Convert to using corners
LE = self.pf.h.gridLeftEdge
RE = self.pf.h.gridRightEdge
vertices = na.array([[LE[:,0],LE[:,1],LE[:,2]],
@@ -508,9 +517,8 @@
self._grids = self.hierarchy.grids[
na.where(na.logical_not(na.all(D<0,axis=0) | na.all(D>0,axis=0) )) ]
+ @cache_mask
def _get_cut_mask(self, grid):
- if self._cut_masks.has_key(grid.id):
- return self._cut_masks[grid.id]
# This is slow. Suggestions for improvement would be great...
ss = grid.ActiveDimensions
D = na.ones(ss) * self._d
@@ -520,8 +528,7 @@
diag_dist = na.sqrt(grid.dx**2.0
+ grid.dy**2.0
+ grid.dz**2.0)
- cm = na.where(na.abs(D) <= 0.5*diag_dist)
- self._cut_masks[grid.id] = cm
+ cm = (na.abs(D) <= 0.5*diag_dist) # Boolean
return cm
def _generate_coords(self):
@@ -558,7 +565,7 @@
def _get_point_indices(self, grid, use_child_mask=True):
k = na.zeros(grid.ActiveDimensions, dtype='bool')
- k[self._get_cut_mask(grid)] = True
+ k = (k | self._get_cut_mask(grid))
if use_child_mask:
pointI = na.where(k & grid.child_mask)
else:
@@ -961,7 +968,7 @@
@restore_grid_state
def _get_point_indices(self, grid, field=None):
k = na.zeros(grid.ActiveDimensions, dtype='bool')
- k[self._get_cut_mask(grid)] = True
+ k = (k | self._get_cut_mask(grid))
pointI = na.where(k & grid.child_mask)
return pointI
@@ -1147,6 +1154,59 @@
def _get_point_indices(self, grid):
return self._indices[grid.id-1]
+class EnzoCylinderBase(Enzo3DData):
+ """
+ We define a disk as have an 'up' vector, a radius and a height.
+ """
+ def __init__(self, center, normal, radius, height, fields=None, pf=None):
+ Enzo3DData.__init__(self, na.array(center), fields, pf)
+ self._norm_vec = na.array(normal)/na.sqrt(na.dot(normal,normal))
+ self.set_field_parameter("height_vector", self._norm_vec)
+ self._height = height
+ self._radius = radius
+ self._d = -1.0 * na.dot(self._norm_vec, self.center)
+ self._cut_masks = {}
+ self._refresh_data()
+
+ def _get_list_of_grids(self):
+ H = na.sum(self._norm_vec.reshape((1,3,1)) * self.pf.h.gridCorners,
+ axis=1) + self._d
+ D = na.sqrt(na.sum((self.pf.h.gridCorners -
+ self.center.reshape((1,3,1)))**2.0,axis=1))
+ R = na.sqrt(D**2.0-H**2.0)
+ self._grids = self.hierarchy.grids[
+ ( (na.any(na.abs(H)<self._height,axis=0))
+ & (na.any(R<self._radius,axis=0)
+ & (na.logical_not((na.all(H>0,axis=0) | (na.all(H<0, axis=0)))) )
+ ) ) ]
+ self._grids = self.hierarchy.grids
+
+ @cache_mask
+ def _get_cut_mask(self, grid):
+ corners = grid._corners.reshape((8,3,1))
+ H = na.sum(self._norm_vec.reshape((1,3,1)) * corners,
+ axis=1) + self._d
+ D = na.sqrt(na.sum((corners -
+ self.center.reshape((1,3,1)))**2.0,axis=1))
+ R = na.sqrt(D**2.0-H**2.0)
+ if na.all(na.abs(H) < self._height, axis=0) \
+ and na.all(R < self._radius, axis=0):
+ cm = na.ones(grid.ActiveDimensions, dtype='bool')
+ else:
+ h = grid['x'] * self._norm_vec[0] \
+ + grid['y'] * self._norm_vec[1] \
+ + grid['z'] * self._norm_vec[2] \
+ + self._d
+ d = na.sqrt(
+ (grid['x'] - self.center[0])**2.0
+ + (grid['y'] - self.center[1])**2.0
+ + (grid['z'] - self.center[2])**2.0
+ )
+ r = na.sqrt(d**2.0-h**2.0)
+ cm = ( (na.abs(h) < self._height)
+ & (r < self._radius))
+ return cm
+
class EnzoRegionBase(Enzo3DData):
"""
EnzoRegions are rectangular prisms of data.
@@ -1164,24 +1224,26 @@
Enzo3DData.__init__(self, center, fields, pf)
self.left_edge = left_edge
self.right_edge = right_edge
+ self._cut_masks = {}
self._refresh_data()
def _get_list_of_grids(self):
self._grids, ind = self.pf.hierarchy.get_box_grids(self.left_edge,
self.right_edge)
+ @cache_mask
def _get_cut_mask(self, grid):
if na.all( (grid._corners < self.right_edge)
& (grid._corners >= self.left_edge)):
- return na.ones(grid.ActiveDimensions, dtype='bool')
- pointI = na.where(\
- ( (grid['x'] < self.right_edge[0])
- & (grid['x'] >= self.left_edge[0])
- & (grid['y'] < self.right_edge[1])
- & (grid['y'] >= self.left_edge[1])
- & (grid['z'] < self.right_edge[2])
- & (grid['z'] >= self.left_edge[2]) ))
- return pointI
+ cm = na.ones(grid.ActiveDimensions, dtype='bool')
+ else:
+ cm = ( (grid['x'] < self.right_edge[0])
+ & (grid['x'] >= self.left_edge[0])
+ & (grid['y'] < self.right_edge[1])
+ & (grid['y'] >= self.left_edge[1])
+ & (grid['z'] < self.right_edge[2])
+ & (grid['z'] >= self.left_edge[2]) )
+ return cm
class EnzoGridCollection(Enzo3DData):
"""
@@ -1198,13 +1260,15 @@
Enzo3DData.__init__(self, center, fields, pf)
self._grids = na.array(grid_list)
self.fields = fields
+ self._cut_masks = {}
self.connection_pool = True
def _get_list_of_grids(self):
pass
+ @cache_mask
def _get_cut_mask(self, grid):
- return na.indices(grid.ActiveDimensions)
+ return na.ones(grid.ActiveDimensions, dtype='bool')
def _get_point_indices(self, grid, use_child_mask=True):
k = na.ones(grid.ActiveDimensions, dtype='bool')
@@ -1239,7 +1303,7 @@
grids.sort(key=lambda x: (x.Level, x.LeftEdge[0], x.LeftEdge[1], x.LeftEdge[2]))
self._grids = na.array(grids)
- @restore_grid_state
+ @restore_grid_state # Pains me not to decorate with cache_mask here
def _get_cut_mask(self, grid, field=None):
# We have the *property* center, which is not necessarily
# the same as the field_parameter
@@ -1250,11 +1314,10 @@
return na.zeros(grid.ActiveDimensions, dtype='bool')
if self._cut_masks.has_key(grid.id):
return self._cut_masks[grid.id]
- pointI = na.where(
- (grid["RadiusCode"]<=self.radius) &
- (grid.child_mask==1) )
- self._cut_masks[grid.id] = pointI
- return pointI
+ cm = ( (grid["RadiusCode"]<=self.radius) &
+ (grid.child_mask==1) )
+ self._cut_masks[grid.id] = cm
+ return cm
class EnzoCoveringGrid(Enzo3DData):
"""
Modified: trunk/yt/lagos/DerivedFields.py
==============================================================================
--- trunk/yt/lagos/DerivedFields.py (original)
+++ trunk/yt/lagos/DerivedFields.py Tue Jan 29 16:38:11 2008
@@ -331,6 +331,8 @@
return data.convert("cm")
add_field("Height", convert_function=_convertHeight,
validators=[ValidateParameter("height_vector")])
+add_field("HeightCode", function=_Height,
+ validators=[ValidateParameter("height_vector")])
def _DynamicalTime(field, data):
"""
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Tue Jan 29 16:38:11 2008
@@ -187,6 +187,7 @@
self.sphere = classobj("EnzoSphere",(EnzoSphereBase,), dd)
self.cutting = classobj("EnzoCuttingPlane",(EnzoCuttingPlaneBase,), dd)
self.ray = classobj("EnzoOrthoRay",(EnzoOrthoRayBase,), dd)
+ self.disk = classobj("EnzoCylinder",(EnzoCylinderBase,), dd)
def __initialize_data_file(self):
"""
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Tue Jan 29 16:38:11 2008
@@ -36,7 +36,6 @@
source.field_parameters = old_params
else:
tr = func(*args, **kwargs)
- #print func.func_name, tr
return tr
return save_state
@@ -80,7 +79,6 @@
used = (used | u)
grid.clear_data()
ub = na.where(used)
- print ub
for field in fields:
if weight:
data[field][ub] /= weight_data[field][ub]
@@ -252,9 +250,6 @@
return
bin_indices_x = na.digitize(sd_x, self[self.x_bin_field])
bin_indices_y = na.digitize(sd_y, self[self.y_bin_field])
- #print bin_indices_x.max(), bin_indices_y.max(),
- #print bin_indices_x.min(), bin_indices_y.min(),
- #print sd_x.size, sd_y.size
# Now we set up our inverse bin indices
return (bin_indices_x, bin_indices_y)
More information about the yt-svn
mailing list