[Yt-svn] yt-commit r391 - trunk/yt/lagos

mturk at wrangler.dreamhost.com mturk at wrangler.dreamhost.com
Wed Apr 2 12:56:02 PDT 2008


Author: mturk
Date: Wed Apr  2 12:56:02 2008
New Revision: 391
URL: http://yt.spacepope.org/changeset/391

Log:
* Changed variable_length to particle_type, as variable_length added unncessary
  ambiguity.  For now, if something is not a grid type, it's a particle.
* 3D objects should now only include particles from within that object.  I did
  a bit of trickery here (and tested it with spheres and regions, but not
  cylinders) and added a mock grid class.  Maybe it could be made faster,
  but this preserved the single cutting mask check.
* Changed a couple of the 3D objects to only check the _is_fully_enclosed
  method and return True if that is True.  This is more DRY-ish.



Modified:
   trunk/yt/lagos/BaseDataTypes.py
   trunk/yt/lagos/BaseGridType.py
   trunk/yt/lagos/DerivedFields.py

Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py	(original)
+++ trunk/yt/lagos/BaseDataTypes.py	Wed Apr  2 12:56:02 2008
@@ -36,12 +36,35 @@
 
 def cache_mask(func):
     def check_cache(self, grid):
-        if not self._cut_masks.has_key(grid.id):
+        if isinstance(grid, FakeGridForParticles):
+            return func(self, grid)
+        elif 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 FakeGridForParticles(object):
+    def __init__(self, grid):
+        self._corners = grid._corners
+        self.field_parameters = {}
+        self.data = {'x':grid['particle_position_x'],
+                     'y':grid['particle_position_y'],
+                     'z':grid['particle_position_z']}
+        self.real_grid = grid
+        self.child_mask = 1
+    def __getitem__(self, field):
+        if field not in self.data.keys():
+            if field == "RadiusCode":
+                center = self.field_parameters['center']
+                tr = na.sqrt( (self['x'] - center[0])**2.0 +
+                              (self['y'] - center[1])**2.0 +
+                              (self['z'] - center[2])**2.0 )
+            else:
+                raise KeyError(field)
+        else: tr = self.data[field]
+        return tr
+
 class EnzoData:
     """
     Generic EnzoData container.  By itself, will attempt to
@@ -443,7 +466,7 @@
         slHERE = tuple(sl)
         sl.reverse()
         slHDF = tuple(sl)
-        if fieldInfo.has_key(field) and fieldInfo[field].variable_length:
+        if fieldInfo.has_key(field) and fieldInfo[field].particle_type:
             return grid[field]
         if not grid.has_key(field):
             conv_factor = 1.0
@@ -547,7 +570,7 @@
         return na.array(coords).swapaxes(0,1)
 
     def _get_data_from_grid(self, grid, field):
-        if not fieldInfo[field].variable_length:
+        if not fieldInfo[field].particle_type:
             pointI = self._get_point_indices(grid)
             if grid[field].size == 1: # dx, dy, dz, cellvolume
                 t = grid[field] * na.ones(grid.ActiveDimensions)
@@ -935,9 +958,10 @@
 
     @restore_grid_state
     def _get_data_from_grid(self, grid, field):
-        if fieldInfo.has_key(field) and fieldInfo[field].variable_length:
+        if field in fieldInfo and fieldInfo[field].particle_type:
+            pointI = self._get_particle_indices(grid)
             try:
-                tr = grid[field]
+                tr = grid[field][pointI].ravel()
             except grid._read_exception:
                 tr = []
             return tr
@@ -996,6 +1020,15 @@
         if use_child_mask: k = (k & grid.child_mask)
         return na.where(k)
 
+    def _get_cut_particle_mask(self, grid):
+        fake_grid = FakeGridForParticles(grid)
+        return self._get_cut_mask(fake_grid)
+
+    def _get_particle_indices(self, grid):
+        k = na.zeros(grid.NumberOfParticles, dtype='bool')
+        k = (k | self._get_cut_particle_mask(grid))
+        return na.where(k)
+
     def extract_region(self, indices):
         """
         Return an ExtractedRegion where the points contained in it are defined
@@ -1149,6 +1182,11 @@
         self._grids = None
         self._refresh_data()
 
+    def _get_cut_particle_mask(self, grid):
+        # Override to provide a warning
+        mylog.warning("Returning all particles from an Extracted Region.  This could be incorrect!")
+        return True
+
     def _get_list_of_grids(self):
         # Okay, so what we're going to want to do is get the pointI from
         # region._get_point_indices(grid) for grid in base_region._grids,
@@ -1222,15 +1260,8 @@
 
     @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')
+        if self._is_fully_enclosed(grid):
+            return True
         else:
             h = grid['x'] * self._norm_vec[0] \
               + grid['y'] * self._norm_vec[1] \
@@ -1276,9 +1307,8 @@
 
     @cache_mask
     def _get_cut_mask(self, grid):
-        if na.all( (grid._corners < self.right_edge)
-                 & (grid._corners >= self.left_edge)):
-            cm = na.ones(grid.ActiveDimensions, dtype='bool')
+        if self._is_fully_enclosed(grid):
+            return True
         else:
             cm = ( (grid['x'] < self.right_edge[0])
                  & (grid['x'] >= self.left_edge[0])
@@ -1355,14 +1385,13 @@
     def _get_cut_mask(self, grid, field=None):
         # We have the *property* center, which is not necessarily
         # the same as the field_parameter
-        corner_radius = na.sqrt(((grid._corners - self.center)**2.0).sum(axis=1))
-        if na.all(corner_radius <= self.radius):
-            return grid.child_mask
-        if self._cut_masks.has_key(grid.id):
+        if self._is_fully_enclosed(grid):
+            return True # We do not want child masking here
+        if not isinstance(grid, FakeGridForParticles) \
+           and grid.id in self._cut_masks:
             return self._cut_masks[grid.id]
-        cm = ( (grid["RadiusCode"]<=self.radius) &
-               (grid.child_mask==1) )
-        self._cut_masks[grid.id] = cm
+        cm = ( (grid["RadiusCode"]<=self.radius) & grid.child_mask )
+        if not isinstance(grid, FakeGridForParticles): self._cut_masks[grid.id] = cm
         return cm
 
 class EnzoCoveringGrid(Enzo3DData):

Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py	(original)
+++ trunk/yt/lagos/BaseGridType.py	Wed Apr  2 12:56:02 2008
@@ -97,7 +97,7 @@
                 try:
                     self[field] = self.readDataFast(field) * conv_factor
                 except self._read_exception:
-                    if field in fieldInfo and fieldInfo[field].variable_length:
+                    if field in fieldInfo and fieldInfo[field].particle_type:
                         self[field] = na.array([],dtype='float64')
                     else: raise
             else:

Modified: trunk/yt/lagos/DerivedFields.py
==============================================================================
--- trunk/yt/lagos/DerivedFields.py	(original)
+++ trunk/yt/lagos/DerivedFields.py	Wed Apr  2 12:56:02 2008
@@ -96,7 +96,7 @@
                  convert_function = None,
                  units = "", projected_units = "",
                  take_log = True, validators = None,
-                 variable_length = False, vector_field=False,
+                 particle_type = False, vector_field=False,
                  line_integral = True,
                  projection_conversion = "cm"):
         self.name = name
@@ -111,7 +111,7 @@
         if not convert_function:
             convert_function = lambda a: 1.0
         self._convert_function = convert_function
-        self.variable_length = variable_length
+        self.particle_type = particle_type
         self.vector_field = vector_field
         self.line_integral = line_integral
         self.projection_conversion = projection_conversion
@@ -297,9 +297,9 @@
     pfunc = particle_func(pf)
     add_field("particle_%s" % pf, function=pfunc,
               validators = [ValidateSpatial(0)],
-              variable_length=True)
+              particle_type=True)
 add_field("particle mass", function=particle_func("particle mass"),
-          validators=[ValidateSpatial(0)], variable_length=True)
+          validators=[ValidateSpatial(0)], particle_type=True)
 
 def _ParticleMass(field, data):
     particles = data["particle mass"] * \
@@ -311,10 +311,10 @@
 def _convertParticleMassMsun(data):
     return data.convert("Density")*((data.convert("cm")**3.0)/1.989e33)
 add_field("ParticleMass", validators=[ValidateSpatial(0)],
-          variable_length=True, convert_function=_convertParticleMass)
+          particle_type=True, convert_function=_convertParticleMass)
 add_field("ParticleMassMsun",
           function=_ParticleMass, validators=[ValidateSpatial(0)],
-          variable_length=True, convert_function=_convertParticleMassMsun)
+          particle_type=True, convert_function=_convertParticleMassMsun)
 
 def _MachNumber(field, data):
     """M{|v|/t_sound}"""



More information about the yt-svn mailing list