[Yt-svn] yt-commit r367 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Jan 24 22:06:52 PST 2008
Author: mturk
Date: Thu Jan 24 22:06:51 2008
New Revision: 367
URL: http://yt.spacepope.org/changeset/367
Log:
Okay, some big changes here.
Fixes and changes:
* Previously, due to an unnoticed bug I introduced, family-trees were not
getting set. I would like to note that this is due to a typo, *not* due to a
problem with the algorithm. :) Mathematical correctness for projections has
been restored, in at least the three datasets I have checked.
* Projections no longer use PointCombine. Hopefully PointCombine.c will not
be long for the SVN-repo.
* Projections have been refactored a bit, and less int-math is used. I think
that the MIPs are not working right now, but I intend to take a look at that
sometime in the near future. Default behavior works.
* ExtractedRegions should be significantly faster now, as I have cached more
than I previously did in terms of indices and so forth. They are still too
slow, so I am going to investigate decorator-based caching to speed things up
even further with less chance of error.
New:
* Added an OrthoRay object, which allows you to cast a ray through the volume.
More on this later, but it is *already* trivial to use this to do proper,
back-to-front line integrals.
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/BaseGridType.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/WeaveStrings.py
trunk/yt/lagos/__init__.py
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Thu Jan 24 22:06:51 2008
@@ -134,6 +134,89 @@
def _generate_field_in_grids(self, fieldName):
pass
+class Enzo1DData(EnzoData):
+ _spatial = False
+ def __init__(self, pf, fields):
+ EnzoData.__init__(self, pf, fields)
+ self._grids = None
+
+ def _generate_field_in_grids(self, field, num_ghost_zones=0):
+ for grid in self._grids:
+ temp = grid[field]
+
+ def _generate_field(self, field):
+ if fieldInfo.has_key(field):
+ # First we check the validator
+ try:
+ fieldInfo[field].check_available(self)
+ except NeedsGridType, ngt_exception:
+ # We leave this to be implementation-specific
+ self._generate_field_in_grids(field, ngt_exception.ghost_zones)
+ return False
+ else:
+ self[field] = fieldInfo[field](self)
+ return True
+ else: # Can't find the field, try as it might
+ raise exceptions.KeyError(field)
+
+
+class EnzoOrthoRayBase(Enzo1DData):
+ def __init__(self, axis, coords, fields=None, pf=None):
+ Enzo1DData.__init__(self, pf, fields)
+ self.axis = axis
+ self.px_ax = x_dict[self.axis]
+ self.py_ax = y_dict[self.axis]
+ self.px_dx = 'd%s'%(axis_names[self.px_ax])
+ self.py_dx = 'd%s'%(axis_names[self.py_ax])
+ self.px, self.py = coords
+ self._refresh_data()
+
+ def _get_list_of_grids(self):
+ y = na.where( (self.px > self.pf.hierarchy.gridLeftEdge[:,self.px_ax])
+ & (self.px < self.pf.hierarchy.gridRightEdge[:,self.px_ax])
+ & (self.py > self.pf.hierarchy.gridLeftEdge[:,self.py_ax])
+ & (self.py < self.pf.hierarchy.gridRightEdge[:,self.py_ax]))
+ self._grids = self.hierarchy.grids[y]
+
+ def get_data(self, fields=None, in_grids=False):
+ if self._grids == None:
+ self._get_list_of_grids()
+ points = []
+ #if not self.has_key('dx'):
+ #self._generate_coords()
+ if not fields:
+ fields_to_get = self.fields
+ else:
+ fields_to_get = ensure_list(fields)
+ mylog.debug("Going to obtain %s (%s)", fields_to_get, self.fields)
+ for field in fields_to_get:
+ if self.data.has_key(field):
+ continue
+ mylog.info("Getting field %s from %s", field, len(self._grids))
+ if field not in self.hierarchy.field_list and not in_grids:
+ if self._generate_field(field):
+ continue # True means we already assigned it
+ self[field] = na.concatenate(
+ [self._get_data_from_grid(grid, field)
+ for grid in self._grids])
+
+ def _get_data_from_grid(self, grid, field):
+ # We are orthogonal, so we can feel free to make assumptions
+ # for the sake of speed.
+ gdx = just_one(grid[self.px_dx])
+ gdy = just_one(grid[self.py_dx])
+ x_coord = int((self.px - grid.LeftEdge[self.px_ax])/gdx)
+ y_coord = int((self.py - grid.LeftEdge[self.py_ax])/gdy)
+ sl = [None,None,None]
+ sl[self.px_ax] = slice(x_coord,x_coord+1,None)
+ sl[self.py_ax] = slice(y_coord,y_coord+1,None)
+ sl[self.axis] = slice(None)
+ if not iterable(grid[field]):
+ gf = grid[field] * na.ones(grid.child_mask[sl].shape)
+ else:
+ gf = grid[field][sl]
+ return gf[na.where(grid.child_mask[sl])]
+
class Enzo2DData(EnzoData):
"""
Class to represent a set of EnzoData that's 2-D in nature, and thus
@@ -473,10 +556,13 @@
def interpolate_discretize(self, *args, **kwargs):
pass
- def _get_point_indices(self, grid):
+ def _get_point_indices(self, grid, use_child_mask=True):
k = na.zeros(grid.ActiveDimensions, dtype='bool')
k[self._get_cut_mask(grid)] = True
- pointI = na.where(k & grid.child_mask)
+ if use_child_mask:
+ pointI = na.where(k & grid.child_mask)
+ else:
+ pointI = na.where(k)
return pointI
class EnzoProjBase(Enzo2DData):
@@ -529,14 +615,14 @@
mylog.debug("Examining level %s", level)
grids = s.levelIndices[level]
num_grids = len(grids)
- RE = s.gridRightEdge[grids].copy()
- LE = s.gridLeftEdge[grids].copy()
+ RE = s.gridRightEdge[grids]
+ LE = s.gridLeftEdge[grids]
for grid in s._grids[grids]:
if (i%1e3) == 0:
mylog.debug("Reading and masking %s / %s", i, len(s._grids))
grid._generate_overlap_masks(self.axis, LE, RE)
grid._overlap_grids[self.axis] = \
- s._grids[grids][na.where(grid.overlap_masks[self.axis] == 1)]
+ s._grids[grids][na.where(grid.overlap_masks[self.axis])]
level_mem[level] += \
grid.ActiveDimensions.prod() / \
grid.ActiveDimensions[self.axis]
@@ -554,7 +640,7 @@
node_name = "%s_%s_%s" % (self.fields[0], self._weight, self.axis)
mylog.info("nodeName: %s", node_name)
projArray = na.array([self['px'], self['py'],
- self['pdx'], self['dy'], self[self.fields[0]]])
+ self['pdx'], self['pdy'], self[self.fields[0]]])
self.hierarchy.save_data(projArray, "/Projections", node_name)
mylog.info("Done serializing...")
@@ -575,23 +661,20 @@
def __project_level(self, level, field):
grids_to_project = self.source.select_grids(level)
zero_out = (level != self._max_level)
- pbar = get_pbar('Projecting level % 2i / % 2i ' \
+ pbar = get_pbar('Projecting level % 2i / % 2i ' \
% (level, self._max_level), len(grids_to_project))
for pi, grid in enumerate(grids_to_project):
grid.retVal = self._project_grid(grid, field, zero_out)
- #grid.retVal = grid._get_projection(self.axis, field, zeroOut,
- #weight=self._weight,
- #func = self.func)
pbar.update(pi)
pbar.finish()
- self.__combine_grids(level) # In-place
+ self.__combine_grids_on_level(level) # In-place
+ if level > 0 and level <= self._max_level:
+ self.__refine_to_level(level) # In-place
all_data = [ [grid.retVal[j] for grid in grids_to_project] for j in range(5)]
for grid in grids_to_project:
cI = na.where(grid.retVal[3]==0) # Where childmask = 0
grid.coarseData = [grid.retVal[j][cI] for j in range(5)]
levelData = [na.concatenate(all_data[i]) for i in range(5)]
- mylog.debug("All done combining and refining with a final %s points",
- levelData[0].shape[0])
dblI = na.where((levelData[0]>-1) & (levelData[3]==1))
if self._weight != None:
weightedData = levelData[2][dblI] / levelData[4][dblI]
@@ -603,38 +686,87 @@
dx = grids_to_project[0].dx * na.ones(len(dblI[0]), dtype='float64')
return [levelData[0][dblI], levelData[1][dblI], weightedData, dx]
- def __combine_grids(self, level):
+ def __setup_weave_dict_combine(self, grid1, grid2):
+ # Recall: x, y, val, mask, weight
+ my_dict = {}
+ my_dict['g1data_x'] = grid1.retVal[0]
+ my_dict['g1data_y'] = grid1.retVal[1]
+ my_dict['g1data_vals'] = grid1.retVal[2]
+ my_dict['g1data_mask'] = grid1.retVal[3]
+ my_dict['g1data_wgt'] = grid1.retVal[4]
+ my_dict['g1points'] = grid1.retVal[0].size
+
+ my_dict['g2data_x'] = grid2.retVal[0]
+ my_dict['g2data_y'] = grid2.retVal[1]
+ my_dict['g2data_vals'] = grid2.retVal[2]
+ my_dict['g2data_mask'] = grid2.retVal[3]
+ my_dict['g2data_wgt'] = grid2.retVal[4]
+ my_dict['g2points'] = grid2.retVal[0].size
+ return my_dict
+
+ def __combine_grids_on_level(self, level):
grids = self.source.select_grids(level)
- pbar = get_pbar('Combining level % 2i / % 2i ' \
+ pbar = get_pbar('Combining level % 2i / % 2i ' \
% (level, self._max_level), len(grids))
# We have an N^2 check, so we try to be as quick as possible
# and to skip as many as possible
for pi, grid1 in enumerate(grids):
pbar.update(pi)
if grid1.retVal[0].shape[0] == 0: continue
- for grid2 in grid1._overlap_grids[self.axis]:
+ for grid2 in grid1._overlap_grids[self.axis].tolist():
if grid2.retVal[0].shape[0] == 0 \
or grid1.id == grid2.id:
continue
- args = grid1.retVal + grid2.retVal + [0]
- PointCombine.CombineData(*args)
+ my_dict = self.__setup_weave_dict_combine(grid1, grid2)
+ weave.inline(ProjectionCombineSameLevel,
+ my_dict.keys(), local_dict=my_dict,
+ compiler='gcc',
+ type_converters=converters.blitz,
+ auto_downcast=0, verbose=2)
goodI = na.where(grid2.retVal[0] > -1)
grid2.retVal = [grid2.retVal[i][goodI] for i in range(5)]
- numRefined = 0
- if level <= self._max_level and level > 0:
- o_grids = grid1.Parent._overlap_grids[self.axis].tolist()
- o_grids.append(grid1.Parent)
- for grid2 in o_grids:
- if grid2.coarseData[0].shape[0] == 0: continue # Already refined
- args = grid1.retVal[:3] + [grid1.retVal[4]] + \
- grid2.coarseData[:3] + [grid2.coarseData[4]] + [2]
- numRefined += PointCombine.RefineCoarseData(*args)
- goodI = na.where(grid2.coarseData[0] > -1)
- grid2.coarseData = [grid2.coarseData[i][goodI] for i in range(5)]
+ pbar.finish()
+
+ def __setup_weave_dict_refine(self, grid1, grid2):
+ # Recall: x, y, val, mask, weight
+ my_dict = {}
+ my_dict['finedata_x'] = grid1.retVal[0]
+ my_dict['finedata_y'] = grid1.retVal[1]
+ my_dict['finedata_vals'] = grid1.retVal[2]
+ my_dict['finedata_wgt'] = grid1.retVal[4]
+ my_dict['fpoints'] = grid1.retVal[0].size
+ my_dict['coarsedata_x'] = grid2.coarseData[0]
+ my_dict['coarsedata_y'] = grid2.coarseData[1]
+ my_dict['coarsedata_vals'] = grid2.coarseData[2]
+ my_dict['coarsedata_wgt'] = grid2.coarseData[4]
+ my_dict['cpoints'] = grid2.coarseData[0].size
+ my_dict['flagged'] = na.zeros(grid2.coarseData[0].size)
+ my_dict['rf'] = int(grid2.dx/grid1.dx)
+ return my_dict
+
+ def __refine_to_level(self, level):
+ grids = self.source.select_grids(level)
+ pbar = get_pbar('Refining to level % 2i / % 2i ' \
+ % (level, self._max_level), len(grids))
+ for pi, grid1 in enumerate(grids):
+ pbar.update(pi)
+ o_grids = grid1.Parent._overlap_grids[self.axis].tolist()
+ for grid2 in o_grids:
+ if grid2.coarseData[0].shape[0] == 0: continue # Already refined
+ my_dict = self.__setup_weave_dict_refine(grid1, grid2)
+ weave.inline(ProjectionRefineCoarseData,
+ my_dict.keys(), local_dict=my_dict,
+ compiler='gcc',
+ type_converters=converters.blitz,
+ auto_downcast=0, verbose=2)
+ goodI = na.where(my_dict['flagged'] == 0)
+ grid2.coarseData = [grid2.coarseData[i][goodI] for i in range(5)]
for grid1 in self.source.select_grids(level-1):
if grid1.coarseData[0].shape[0] != 0:
- mylog.error("Something fucked up, and %s still has %s points of data",
+ mylog.error("Something messed up, and %s still has %s points of data",
grid1, grid1.coarseData[0].size)
+ print grid1.coarseData[0]
+ raise ValueError(grid1)
pbar.finish()
@time_execution
@@ -651,7 +783,7 @@
self['py'] = (all_data[1,:]+0.5) * self['pdx']
self['pdx'] *= 0.5
self['pdy'] = self['pdx'].copy()
- self.data[field] = all_data[2,:]
+ self[field] = all_data[2,:]
# Now, we should clean up after ourselves...
[grid.clear_data() for grid in s._grids]
@@ -660,45 +792,41 @@
masked_data = self._get_data_from_grid(grid, field)
weight_data = na.ones(masked_data.shape)
else:
- masked_data = self._get_data_from_grid(grid, field) \
- * self._get_data_from_grid(grid, self._weight)
weight_data = self._get_data_from_grid(grid, self._weight)
- if zero_out:
+ masked_data = self._get_data_from_grid(grid, field) * weight_data
+ if zero_out: # Not sure how to do this here
masked_data[grid.child_indices] = 0
weight_data[grid.child_indices] = 0
- dl = grid['d%s' % axis_names[self.axis]]
- dx = grid['d%s' % axis_names[x_dict[self.axis]]]
- dy = grid['d%s' % axis_names[y_dict[self.axis]]]
- full_proj = self.func(masked_data,self.axis)*dl
- weight_proj = self.func(weight_data,self.axis)*dl
- used_data = self._get_cut_mask(grid)
- used_points = na.where(self.func(used_data, self.axis) > 0)
+ #@todo: Fix the func to set up a path length too
+ dl = just_one(grid['d%s' % axis_names[self.axis]])
+ dx = just_one(grid['d%s' % axis_names[x_dict[self.axis]]])
+ dy = just_one(grid['d%s' % axis_names[y_dict[self.axis]]])
+ full_proj = self.func(masked_data,axis=self.axis)*dl
+ weight_proj = self.func(weight_data,axis=self.axis)*dl
+ used_data = self._get_points_in_region(grid)
+ used_points = na.where(na.logical_or.reduce(used_data, self.axis))
if zero_out:
subgrid_mask = na.logical_and.reduce(grid.child_mask, self.axis).astype('int64')
else:
subgrid_mask = na.ones(full_proj.shape, dtype='int64')
xind, yind = [arr[used_points].ravel() for arr in na.indices(full_proj.shape)]
- xpoints = xind + na.rint(grid.LeftEdge[x_dict[self.axis]]/dx).astype('int64')
- ypoints = yind + na.rint(grid.LeftEdge[y_dict[self.axis]]/dy).astype('int64')
+ start_index = grid.get_global_startindex()
+ xpoints = xind + (start_index[x_dict[self.axis]])
+ ypoints = yind + (start_index[y_dict[self.axis]])
return [xpoints, ypoints,
full_proj[used_points].ravel(),
subgrid_mask[used_points].ravel(),
weight_proj[used_points].ravel()]
- def _get_cut_mask(self, grid):
- tr = na.zeros(grid.ActiveDimensions, dtype='bool')
- tr[self.source._get_cut_mask(grid)] = True
- return tr
-
- def _get_used_points(self, grid):
- pointI = self.source._get_point_indices(grid)
- bad_points = na.zeros(grid.ActiveDimensions)
- bad_points[pointI] = 1.0
- return bad_points
+ def _get_points_in_region(self, grid):
+ pointI = self.source._get_point_indices(grid, use_child_mask=False)
+ point_mask = na.zeros(grid.ActiveDimensions)
+ point_mask[pointI] = 1.0
+ return point_mask
@restore_grid_state
def _get_data_from_grid(self, grid, field):
- bad_points = self._get_used_points(grid)
+ bad_points = self._get_points_in_region(grid)
d = grid[field] * bad_points
return d
@@ -726,6 +854,7 @@
self.set_field_parameter("center",center)
self.coords = None
self.dx = None
+ self._grids = None
def _generate_coords(self):
mylog.info("Generating coords for %s grids", len(self._grids))
@@ -757,7 +886,8 @@
return tr
def get_data(self, fields=None, in_grids=False):
- self._get_list_of_grids()
+ if self._grids == None:
+ self._get_list_of_grids()
points = []
#if not self.has_key('dx'):
#self._generate_coords()
@@ -923,7 +1053,7 @@
__del_gridLevels)
def extract_connected_sets(self, field, num_levels, min_val, max_val,
- log_space=True):
+ log_space=True, cumulative=True):
"""
This function will create a set of contour objects, defined
by having connected cell structures, which can then be
@@ -938,7 +1068,11 @@
contours = {}
for level in range(num_levels):
contours[level] = {}
- cids = identify_contours(self, field, cons[level], cons[level+1])
+ if cumulative:
+ mv = max_val
+ else:
+ mv = cons[level+1]
+ cids = identify_contours(self, field, cons[level], mv)
for cid, cid_ind in cids.items():
contours[level][cid] = self.extract_region(cid_ind)
return cons, contours
@@ -980,12 +1114,14 @@
fields=None, pf=base_region.pf)
self._base_region = base_region
self._base_indices = indices
+ self._grids = None
self._refresh_data()
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,
# and then construct an array of those, which we will select along indices.
+ if self._grids != None: return
grid_vals, xi, yi, zi = [], [], [], []
for grid in self._base_region._grids:
xit,yit,zit = self._base_region._get_point_indices(grid)
@@ -1008,9 +1144,8 @@
zi[self._base_indices][ind_ind]])
self._grids = self._base_region.pf.h.grids[self._indices.keys()]
- def _get_cut_mask(self, grid):
- x,y,z = self._indices[grid.id-1]
- return (x,y,z)
+ def _get_point_indices(self, grid):
+ return self._indices[grid.id-1]
class EnzoRegionBase(Enzo3DData):
"""
@@ -1039,13 +1174,13 @@
if na.all( (grid._corners < self.right_edge)
& (grid._corners >= self.left_edge)):
return na.ones(grid.ActiveDimensions, dtype='bool')
- pointI = \
+ 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]) )
+ & (grid['z'] >= self.left_edge[2]) ))
return pointI
class EnzoGridCollection(Enzo3DData):
@@ -1069,7 +1204,14 @@
pass
def _get_cut_mask(self, grid):
- return na.ones(grid.ActiveDimensions, dtype='bool')
+ return na.indices(grid.ActiveDimensions)
+
+ def _get_point_indices(self, grid, use_child_mask=True):
+ k = na.ones(grid.ActiveDimensions, dtype='bool')
+ if use_child_mask:
+ k[grid.child_indices] = False
+ pointI = na.where(k == True)
+ return pointI
class EnzoSphereBase(Enzo3DData):
"""
@@ -1108,8 +1250,9 @@
return na.zeros(grid.ActiveDimensions, dtype='bool')
if self._cut_masks.has_key(grid.id):
return self._cut_masks[grid.id]
- pointI = ((grid["RadiusCode"]<=self.radius) &
- (grid.child_mask==1))
+ pointI = na.where(
+ (grid["RadiusCode"]<=self.radius) &
+ (grid.child_mask==1) )
self._cut_masks[grid.id] = pointI
return pointI
Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py (original)
+++ trunk/yt/lagos/BaseGridType.py Thu Jan 24 22:06:51 2008
@@ -50,6 +50,7 @@
self.data = {}
self.field_parameters = {}
self.fields = []
+ self.start_index = None
self.id = id
if hierarchy: self.hierarchy = hierarchy
if filename: self.set_filename(filename)
@@ -156,11 +157,40 @@
self.dy = self.Parent.dy/2.0
self.dz = self.Parent.dz/2.0
ParentLeftIndex = na.rint((self.LeftEdge-self.Parent.LeftEdge)/self.Parent.dx)
+ self.start_index = 2*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
self.LeftEdge = self.Parent.LeftEdge + self.Parent.dx * ParentLeftIndex
- self.RightEdge = self.LeftEdge + self.ActiveDimensions*self.dx
+ self.RightEdge = self.LeftEdge + \
+ self.ActiveDimensions*na.array([self.dx,self.dy,self.dz])
self.hierarchy.gridDxs[self.id-1,0] = self.dx
self.hierarchy.gridDys[self.id-1,0] = self.dy
self.hierarchy.gridDzs[self.id-1,0] = self.dz
+ self.hierarchy.gridLeftEdge[self.id-1,:] = self.LeftEdge
+ self.hierarchy.gridRightEdge[self.id-1,:] = self.RightEdge
+ self.hierarchy.gridCorners[:,:,self.id-1] = na.array([ # Unroll!
+ [self.LeftEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.RightEdge[1], self.LeftEdge[2]],
+ [self.RightEdge[0], self.RightEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.RightEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.LeftEdge[1], self.RightEdge[2]],
+ [self.RightEdge[0], self.LeftEdge[1], self.RightEdge[2]],
+ [self.LeftEdge[0], self.RightEdge[1], self.LeftEdge[2]],
+ ], dtype='float64')
+ self.__child_mask = None
+ self.__child_indices = None
+ self._setup_dx()
+
+ def get_global_startindex(self):
+ if self.start_index != None:
+ return self.start_index
+ if self.Parent == None:
+ start_index = self.LeftEdge / na.array([self.dx, self.dy, self.dz])
+ return na.rint(start_index).astype('int64').ravel()
+ pdx = na.array([self.Parent.dx, self.Parent.dy, self.Parent.dz]).ravel()
+ start_index = (self.Parent.get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
+ self.start_index = (start_index*2).astype('int64').ravel()
+ return self.start_index
def _generate_overlap_masks(self, axis, LE, RE):
"""
@@ -183,7 +213,7 @@
cond3 = self.RightEdge[y] >= LE[:,y]
cond4 = self.LeftEdge[y] <= RE[:,y]
self.overlap_masks[axis]=na.logical_and(na.logical_and(cond1, cond2), \
- na.logical_and(cond3, cond4))
+ na.logical_and(cond3, cond4))
def __repr__(self):
return "Grid_%04i" % (self.id)
@@ -217,7 +247,7 @@
@param field: field to check
@type field: string
"""
- coord=nd.maximum_position(self[field])
+ coord=nd.maximum_position(self[field]*self.child_mask)
val = self[field][coord]
return val, coord
@@ -393,13 +423,11 @@
# Now let's get our overlap
si = [None]*3
ei = [None]*3
- startIndex = ((child.LeftEdge - self.LeftEdge)/self.dx)
- endIndex = ((child.RightEdge - self.LeftEdge)/self.dx)
+ startIndex = na.rint((child.LeftEdge - self.LeftEdge)/self.dx)
+ endIndex = na.rint((child.RightEdge - self.LeftEdge)/self.dx)
for i in xrange(3):
- si[i] = na.rint(startIndex[i])
- ei[i] = na.rint(endIndex[i])
- si[i] = si[i]
- ei[i] = ei[i]
+ si[i] = startIndex[i]
+ ei[i] = endIndex[i]
self.__child_mask[si[0]:ei[0], si[1]:ei[1], si[2]:ei[2]] = 0
self.__child_indices = na.where(self.__child_mask==0)
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Thu Jan 24 22:06:51 2008
@@ -186,6 +186,7 @@
self.covering_grid = classobj("EnzoCoveringGrid",(EnzoCoveringGrid,), dd)
self.sphere = classobj("EnzoSphere",(EnzoSphereBase,), dd)
self.cutting = classobj("EnzoCuttingPlane",(EnzoCuttingPlaneBase,), dd)
+ self.ray = classobj("EnzoOrthoRay",(EnzoOrthoRayBase,), dd)
def __initialize_data_file(self):
"""
@@ -367,9 +368,9 @@
self.gridLevels[secondGrid] = self.gridLevels[firstGrid] + 1
elif m.group(2) == "This":
parent = self.gridReverseTree[firstGrid]
- if parent:
+ if parent and parent > -1:
self.gridTree[parent-1].append(self.grids[secondGrid])
- self.gridReverseTree[secondGrid] = parent
+ self.gridReverseTree[secondGrid] = parent
self.grids[secondGrid].Level = self.grids[firstGrid].Level
self.gridLevels[secondGrid] = self.gridLevels[firstGrid]
pTree = [ [ grid.id - 1 for grid in self.gridTree[i] ] for i in range(self.num_grids) ]
@@ -471,7 +472,8 @@
mylog.debug("Done flushing to grids")
if ytcfg.getboolean("lagos","ReconstructHierarchy") == True:
mylog.debug("Reconstructing hierarchy")
- for grid in self.grids:
+ for level in range(self.maxLevel+1):
+ for grid in self.select_grids(level):
if grid.Parent: grid._guess_properties_from_parent()
def __select_level(self, level):
@@ -635,8 +637,8 @@
pos[0] += 0.5*maxGrid.dx
pos[1] += 0.5*maxGrid.dx
pos[2] += 0.5*maxGrid.dx
- mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s", \
- maxVal, pos[0], pos[1], pos[2], maxGrid, maxGrid.Level)
+ mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
+ maxVal, pos[0], pos[1], pos[2], maxGrid, maxGrid.Level, mc)
self.center = pos
# This probably won't work for anyone else
self.bulkVelocity = (maxGrid["x-velocity"][maxCoord], \
Modified: trunk/yt/lagos/WeaveStrings.py
==============================================================================
--- trunk/yt/lagos/WeaveStrings.py (original)
+++ trunk/yt/lagos/WeaveStrings.py Thu Jan 24 22:06:51 2008
@@ -71,52 +71,47 @@
}
"""
-ProjectionRefineCoarseData = \
-"""
-/*
-fpoints, finedata_x, finedata_y, finedata_vals, finedata_wgt,
-cpoints, coarsedata_x, coarsedata_y, coarsedata_vals, coarsedata_wgt,
-refinementFactor, totalRefined
-
-#define MIPCOMB(A,B) ((A) > (B) ? (A) : (B))
-#define SUMCOMB(A,B) (A + B)
-
-long fi, ci;
-
-int flagged[cpoints], tr = 0;
-
-long double rf = refinementFactor;
+ProjectionCombineSameLevel = \
+r"""
+ int g1i, g2i;
-npy_int64 coarseCell_x, coarseCell_y;
+ for (g1i = 0; g1i < g1points; g1i++) {
+ if (g1data_x(g1i) < 0) continue;
+ for (g2i = 0; g2i < g2points; g2i++) {
+ if (g2data_x(g2i) < 0) continue;
+ if ((g1data_x(g1i) == g2data_x(g2i)) &&
+ (g1data_y(g1i) == g2data_y(g2i))) {
+ g1data_vals(g1i) += g2data_vals(g2i);
+ g1data_wgt(g1i) += g2data_wgt(g2i);
+ g1data_mask(g1i) = ((g1data_mask(g1i)) && (g2data_mask(g2i)));
+ g2data_x(g2i) = -1;
+ break; // We map to one g1data point AT MOST
+ }
+ }
+ }
+"""
-for (ci=0; ci<cpoints; ci++) flagged[ci]=0;
+ProjectionRefineCoarseData = \
+r"""
+ int ci, fi, x_off, y_off;
+ long int fine_x, fine_y;
-for (fi = 0; fi < fpoints; fi++) {
- coarseCell_x = floorl((finedata_x(fi))/rf);
- coarseCell_y = floorl((finedata_y(fi))/rf);
for (ci = 0; ci < cpoints; ci++) {
- if ((coarseCell_x == coarsedata_x(ci)) &&
- (coarseCell_y == coarsedata_y(ci))) {
- tr += 1;
- finedata_vals(fi) = %(COMBTYPE)sCOMB(coarsedata_vals[ci], finedata_vals[fi]);
- finedata_wgt(fi) = %(COMBTYPE)sCOMB(coarsedata_wgt[ci], finedata_wgt[fi]);
- flagged[ci] = 1;
- break; // Each fine cell maps to one and only one coarse
- // cell
+ for (x_off = 0; x_off < rf; x_off++) {
+ for (y_off = 0; y_off < rf; y_off++) {
+ fine_x = 2*coarsedata_x(ci) + x_off;
+ fine_y = 2*coarsedata_y(ci) + y_off;
+ for (fi = 0; fi < fpoints; fi++) {
+ if ((fine_x == finedata_x(fi)) &&
+ (fine_y == finedata_y(fi))) {
+ finedata_vals(fi) += coarsedata_vals(ci);
+ finedata_wgt(fi) += coarsedata_wgt(ci);
+ flagged(ci) += 1;
+ }
+ }
}
+ }
}
-}
-*totalRefined = tr;
-
-for (ci=0; ci<cpoints; ci++) {
- if (flagged[ci]==1) {
- coarsedata_x(ci)=-1;
- coarsedata_y(ci)=-1;
- coarsedata_vals(ci)=0.0;
- coarsedata_wgt(ci)=0.0;
- }
-}
-*/
"""
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Thu Jan 24 22:06:51 2008
@@ -59,7 +59,6 @@
if ytcfg.getboolean("lagos","usefortran"):
try:
import EnzoFortranRoutines
- import EnzoFortranWrapper
except ImportError:
pass
@@ -72,6 +71,7 @@
from DerivedQuantities import quantityInfo
from DataReadingFuncs import *
from ClusterFiles import *
+from ContourFinder import *
from BaseDataTypes import *
from BaseGridType import *
from EnzoRateData import *
@@ -90,4 +90,4 @@
# We assume that it is with respect to the $HOME/.yt directory
execfile(os.path.expanduser("~/.yt/%s" % my_plugin_name))
-log_fields = [] # @todo: GET RID OF THIS
+log_fields = [] # @todo: GET RID OF THIS
\ No newline at end of file
More information about the yt-svn
mailing list