# HG changeset patch # User Matthew Turk # Date 1247616167 14400 # Branch raytrace # Node ID 04e34eae18f4b98d5e67ceddaae1b8abef194990 # Parent 4f34f2d7c72231a67bce34389f475939d89ca605 raytest.py passes now with the epsilon added to the intersection point. diff -r 4f34f2d7c72231a67bce34389f475939d89ca605 -r 04e34eae18f4b98d5e67ceddaae1b8abef194990 yt/lagos/BaseDataTypes.py --- a/yt/lagos/BaseDataTypes.py Tue Jul 14 14:42:53 2009 -0600 +++ b/yt/lagos/BaseDataTypes.py Tue Jul 14 20:02:47 2009 -0400 @@ -458,16 +458,17 @@ i1 = (i+1) % 3 i2 = (i+2) % 3 vs = self._get_line_at_coord(LE[:,i], i) - p = p | ( ( (LE[:,i1] < vs[:,i1]) & (RE[:,i1] > vs[:,i1]) ) \ - & ( (LE[:,i2] < vs[:,i2]) & (RE[:,i2] > vs[:,i2]) ) ) + p = p | ( ( (LE[:,i1] <= vs[:,i1]) & (RE[:,i1] >= vs[:,i1]) ) \ + & ( (LE[:,i2] <= vs[:,i2]) & (RE[:,i2] >= vs[:,i2]) ) ) vs = self._get_line_at_coord(RE[:,i], i) - p = p | ( ( (LE[:,i1] < vs[:,i1]) & (RE[:,i1] > vs[:,i1]) ) \ - & ( (LE[:,i2] < vs[:,i2]) & (RE[:,i2] > vs[:,i2]) ) ) - p = p | ( na.all( LE < self.start_point, axis=1 ) - & na.all( RE > self.start_point, axis=1 ) ) - p = p | ( na.all( LE < self.end_point, axis=1 ) - & na.all( RE > self.end_point, axis=1 ) ) - self._grids = self.hierarchy.grids.copy()#[p] + p = p | ( ( (LE[:,i1] <= vs[:,i1]) & (RE[:,i1] >= vs[:,i1]) ) \ + & ( (LE[:,i2] <= vs[:,i2]) & (RE[:,i2] >= vs[:,i2]) ) ) + p = p | ( na.all( LE <= self.start_point, axis=1 ) + & na.all( RE >= self.start_point, axis=1 ) ) + p = p | ( na.all( LE <= self.end_point, axis=1 ) + & na.all( RE >= self.end_point, axis=1 ) ) + #self._grids = self.hierarchy.grids[p] + self._grids = self.hierarchy.grids.copy() def _get_line_at_coord(self, v, index): # t*self.vec + self.start_point = self.end_point diff -r 4f34f2d7c72231a67bce34389f475939d89ca605 -r 04e34eae18f4b98d5e67ceddaae1b8abef194990 yt/lagos/RTIntegrator.pyx --- a/yt/lagos/RTIntegrator.pyx Tue Jul 14 14:42:53 2009 -0600 +++ b/yt/lagos/RTIntegrator.pyx Tue Jul 14 20:02:47 2009 -0400 @@ -103,69 +103,63 @@ tr = (right_edge[i] - u[i])/v[i] if (left_edge[x] <= (u[x] + tl*v[x]) <= right_edge[x]) and \ (left_edge[y] <= (u[y] + tl*v[y]) <= right_edge[y]) and \ - (0 <= tl < intersect_t): + (0.0 <= tl < intersect_t): intersect_t = tl if (left_edge[x] <= (u[x] + tr*v[x]) <= right_edge[x]) and \ (left_edge[y] <= (u[y] + tr*v[y]) <= right_edge[y]) and \ - (0 <= tr < intersect_t): + (0.0 <= tr < intersect_t): intersect_t = tr # if fully enclosed if (left_edge[0] <= u[0] <= right_edge[0]) and \ (left_edge[1] <= u[1] <= right_edge[1]) and \ (left_edge[2] <= u[2] <= right_edge[2]): - intersect_t = 0 + intersect_t = 0.0 if not (0 <= intersect_t <= 1): return # Now get the indices of the intersection intersect = u + intersect_t * v + cdef int ncells = 0 for i in range(3): - cur_ind[i] = np.floor((intersect[i] - left_edge[i])/dx[i]) - if cur_ind[i] == grid_mask.shape[i]: cur_ind[i] = grid_mask.shape[i] - 1 + cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i]) tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i] + if cur_ind[i] == grid_mask.shape[i] and step[i] < 0: + cur_ind[i] = grid_mask.shape[i] - 1 if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i] if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i] - tdelta[i] = step[i]*(dx[i]/v[i]) + tdelta[i] = np.abs((dx[i]/v[i])) # The variable intersect contains the point we first pierce the grid enter_t = intersect_t while 1: - if not (0 <= cur_ind[0] < grid_mask.shape[0]) or \ - not (0 <= cur_ind[1] < grid_mask.shape[1]) or \ - not (0 <= cur_ind[2] < grid_mask.shape[2]): + if (not (0 <= cur_ind[0] < grid_mask.shape[0])) or \ + (not (0 <= cur_ind[1] < grid_mask.shape[1])) or \ + (not (0 <= cur_ind[2] < grid_mask.shape[2])): break # Note that we are calculating t on the fly, but we get *negative* t # values from what they should be. # If we've reached t = 1, we are done. - if (tmax[0] >= 1.0) and (tmax[1] >= 1.0) and (tmax[2] >= 1.0): + grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 + if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0): grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = 1.0 - enter_t - if grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] > dt_tolerance: - grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 break + ncells += 1 if tmax[0] < tmax[1]: if tmax[0] < tmax[2]: grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[0] - enter_t - if grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] > dx[0] * dt_tolerance: - grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 enter_t = tmax[0] tmax[0] += tdelta[0] cur_ind[0] += step[0] else: grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t - if grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] > dx[2] * dt_tolerance: - grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 enter_t = tmax[2] tmax[2] += tdelta[2] cur_ind[2] += step[2] else: if tmax[1] < tmax[2]: grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[1] - enter_t - if grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] > dx[1] * dt_tolerance: - grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 enter_t = tmax[1] tmax[1] += tdelta[1] cur_ind[1] += step[1] else: grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t - if grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] > dx[2] * dt_tolerance: - grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1 enter_t = tmax[2] tmax[2] += tdelta[2] cur_ind[2] += step[2]