[yt-svn] commit/yt-3.0: 10 new changesets
Bitbucket
commits-noreply at bitbucket.org
Tue Sep 25 03:20:45 PDT 2012
10 new commits in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/changeset/8866b8b706e6/
changeset: 8866b8b706e6
branch: yt-3.0
user: scopatz
date: 2012-09-05 01:14:00
summary: Starting an alternative ray tracer module.
affected #: 2 files
diff -r c0f5e4d39a1ad71e2c4f97a70dda62431fd3730f -r 8866b8b706e6e9ec21a456373e1b13c16345077f yt/utilities/lib/alt_ray_tracers.pyx
--- /dev/null
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -0,0 +1,190 @@
+"""Ray tracing algorithms for rays in non-cartesian coordinate systems.
+
+Author: Anthony Scopatz <scopatz at gmail.com>
+Affiliation: University of Chicago
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Anthony Scopatz. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+cimport libc.math as math
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def _cyl2cart(np.ndarray[np.float64_t, ndim=2] x):
+ """Converts points in cylindrical coordinates to cartesian points."""
+ # NOTE this should be removed once the coord interface comes online
+ cdef int i, I
+ cdef np.ndarray[np.float64_t, ndim=2] xcart
+ I = x.shape[0]
+ xcart = np.empty((I, x.shape[1]), dtype='float64')
+ for i in range(I):
+ xcart[i,0] = x[i,0] * math.cos(x[i,2])
+ xcart[i,1] = x[i,0] * math.sin(x[i,2])
+ xcart[i,2] = x[i,1]
+ return xcart
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def _cart_intersect(np.ndarray[np.float64_t, ndim=1] a,
+ np.ndarray[np.float64_t, ndim=1] b,
+ np.ndarray[np.float64_t, ndim=2] c,
+ np.ndarray[np.float64_t, ndim=2] d):
+ """Finds the times and locations of the lines defined by a->b and
+ c->d in cartesian space. a and b must be 1d points. c and d must
+ be 2d arrays of equal shape whose second dim is the same length as
+ a and b.
+ """
+ cdef int i, I
+ cdef np.ndarray[np.float64_t, ndim=1] r, s, t
+ cdef np.ndarray[np.float64_t, ndim=2] loc
+ I = c.shape[0]
+ shape = (I, c.shape[1])
+ t = np.empty(I, dtype='float64')
+ loc = np.empty(shape, dtype='float64')
+
+ r = b - a
+ for i in range(I):
+ s = d[i] - c[i]
+ t[i] = (np.cross(c[i] - a, s).sum()) / (np.cross(r, s).sum())
+ loc[i] = a + t[i]*r
+ return t, loc
+
+
+"""\
+D = F - E
+
+Ecart = na.array((E[0]*na.cos(E[2]), E[0]*na.sin(E[2]), E[1]))
+Fcart = na.array((F[0]*na.cos(F[2]), F[0]*na.sin(F[2]), F[1]))
+Dcart = Fcart - Ecart
+
+# <codecell>
+
+rleft = pf.h.grid_left_edge[:,0]
+rright = pf.h.grid_right_edge[:,0]
+zleft = pf.h.grid_left_edge[:,1]
+zright = pf.h.grid_right_edge[:,1]
+
+a = (Dcart**2)[:2].sum()
+b = (2*Dcart*Ecart)[:2].sum()
+cleft = (Ecart**2)[:2].sum() - rleft**2
+cright = (Ecart**2)[:2].sum() - rright**2
+
+tmleft = (-b - na.sqrt(b**2 - 4*a*cleft)) / (2*a)
+tpleft = (-b + na.sqrt(b**2 - 4*a*cleft)) / (2*a)
+tmright = (-b - na.sqrt(b**2 - 4*a*cright)) / (2*a)
+tpright = (-b + na.sqrt(b**2 - 4*a*cright)) / (2*a)
+
+tmmright = np.logical_and(~np.isnan(tmright), rright <= E[0])
+tpmright = np.logical_and(~np.isnan(tpright), rright <= F[0])
+
+tmmleft = np.logical_and(~np.isnan(tmleft), rleft <= E[0])
+tpmleft = np.logical_and(~np.isnan(tpleft), rleft <= F[0])
+
+# <codecell>
+
+ind = np.unique(np.concatenate([np.argwhere(tmmleft).flat, np.argwhere(tpmleft).flat, np.argwhere(tmmright).flat, np.argwhere(tpmright).flat,]))
+
+thetaleft = np.arctan2((Ecart[1] + tmleft[ind]*Dcart[1]), (Ecart[0] + tmleft[ind]*Dcart[0]))
+nans = np.isnan(thetaleft)
+thetaleft[nans] = np.arctan2((Ecart[1] + tpleft[ind[nans]]*Dcart[1]), (Ecart[0] + tpleft[ind[nans]]*Dcart[0]))
+
+thetaright = np.arctan2((Ecart[1] + tmright[ind]*Dcart[1]), (Ecart[0] + tmright[ind]*Dcart[0]))
+nans = np.isnan(thetaleft)
+thetaright[nans] = np.arctan2((Ecart[1] + tpright[ind[nans]]*Dcart[1]), (Ecart[0] + tpright[ind[nans]]*Dcart[0]))
+
+#thetaleft += np.pi*3/2
+#thetaright += np.pi*3/2
+
+if 0 == len(ind):
+ print "Ind len zero"
+ I = len(rleft)
+ ind = np.arange(I)
+ thetaleft = np.empty(I)
+ thetaleft.fill(E[2])
+ thetaright = np.empty(I)
+ thetaleft.fill(F[2])
+
+# <codecell>
+
+a = E
+b = F
+
+#c = np.array([rleft[ind], zleft[ind], thetaleft]).T
+#d = np.array([rleft[ind], zright[ind], thetaleft]).T
+
+#c = np.array([rright[ind], zleft[ind], thetaright]).T
+#d = np.array([rright[ind], zright[ind], thetaright]).T
+
+c = np.array([rleft[ind], zright[ind], thetaleft]).T
+c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
+c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
+c = np.append(c, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
+c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
+c = np.append(c, np.array([rright[ind], zright[ind], thetaleft]).T, axis=0)
+c = np.append(c, np.array([rleft[ind], zright[ind], thetaleft]).T, axis=0)
+c = np.append(c, np.array([rright[ind], zleft[ind], thetaleft]).T, axis=0)
+
+d = np.array([rright[ind], zright[ind], thetaright]).T
+d = np.append(d, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
+d = np.append(d, np.array([rleft[ind], zright[ind], thetaleft]).T, axis=0)
+d = np.append(d, np.array([rright[ind], zright[ind], thetaright]).T, axis=0)
+d = np.append(d, np.array([rleft[ind], zleft[ind], thetaright]).T, axis=0)
+d = np.append(d, np.array([rright[ind], zright[ind], thetaright]).T, axis=0)
+d = np.append(d, np.array([rleft[ind], zright[ind], thetaright]).T, axis=0)
+d = np.append(d, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
+
+origind = ind
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+ind = np.append(ind, origind)
+
+tsec, intsec = intersect(a, b, c, d)
+print tsec
+print np.isnan(tsec).all(), len(tsec)
+tmask = np.logical_and(0.0<=tsec, tsec<=1.0)
+tsec, utmask = np.unique(tsec[tmask], return_index=True)
+print tsec.max(), tsec.ptp(), len(utmask)
+ind = ind[tmask][utmask]
+xyz = intsec[tmask][utmask]
+s = np.sqrt(((xyz - Ecart)**2).sum(axis=1))
+s, smask = np.unique(s, return_index=True)
+ind = ind[smask]
+xyz = xyz[smask]
+t = s/np.sqrt((Dcart**2).sum())
+si = s.argsort()
+s = s[si]
+t = t[si]
+ind = ind[si]
+xyz = xyz[si]
+rztheta = np.array([np.sqrt(xyz[:,0]**2 + xyz[:,1]**2), xyz[:,2], np.arctan2(xyz[:,1], xyz[:,0])]).T
+rztheta[:,2] = 0.0 + (rztheta[:,2] - np.pi*3/2)%(2*np.pi)
+
+
+
+
+"""
diff -r c0f5e4d39a1ad71e2c4f97a70dda62431fd3730f -r 8866b8b706e6e9ec21a456373e1b13c16345077f yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -135,6 +135,9 @@
config.add_extension("Interpolators",
["yt/utilities/lib/Interpolators.pyx"],
libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ config.add_extension("alt_ray_tracers",
+ ["yt/utilities/lib/alt_ray_tracers.pyx"],
+ libraries=["m"], depends=[])
config.add_extension("marching_cubes",
["yt/utilities/lib/marching_cubes.pyx",
"yt/utilities/lib/FixedInterpolator.c"],
https://bitbucket.org/yt_analysis/yt-3.0/changeset/2020eb46759b/
changeset: 2020eb46759b
branch: yt-3.0
user: scopatz
date: 2012-09-05 21:22:54
summary: Partway through ray tracer.
affected #: 1 file
diff -r 8866b8b706e6e9ec21a456373e1b13c16345077f -r 2020eb46759bd8aa14da171111edad63c6aa1469 yt/utilities/lib/alt_ray_tracers.pyx
--- a/yt/utilities/lib/alt_ray_tracers.pyx
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -71,59 +71,117 @@
return t, loc
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def clyindrical_ray_trace(np.ndarray[np.float64_t, ndim=1] p1,
+ np.ndarray[np.float64_t, ndim=1] p2,
+ np.ndarray[np.float64_t, ndim=2] left_edges,
+ np.ndarray[np.float64_t, ndim=2] right_edges)
+
+ """Computes straight (cartesian) rays being traced through a
+ cylindrical geometry.
+
+ Parameters
+ ----------
+ p1 : length 3 float ndarray
+ start point for ray
+ p2 : length 3 float ndarray
+ stop point for ray
+ left_edges : 2d ndarray
+ left edges of grid cells
+ right_edges : 2d ndarray
+ right edges of grid cells
+
+ Returns
+ -------
+ t : 1d float ndarray
+ ray parametric time on range [0,1]
+ s : 1d float ndarray
+ ray parametric distance on range [0,len(ray)]
+ rztheta : 2d float ndarray
+ ray grid cell intersections in cylidrical coordinates
+ inds : 1d int ndarray
+ indexes into the grid cells which the ray crosses in order.
+
+ """
+ int i, I
+ cdef np.float64_t a, b,b2, twoa
+ cdef np.ndarray[np.float64_t, ndim=1] dp, p1cart, p2cart, dpcart, t, s,
+ rleft, rright, zleft, zright,
+ cleft, cright, thetaleft, thetaright,
+ tmleft, tpleft, tmright, tpright
+ cdef np.ndarray[np.int32_t, ndim=1] inds
+ cdef np.ndarray[np.float64_t, ndim=2] rztheta, ptemp
+
+ # set up points
+ dp = p2 - p1
+ ptemp = np.array([p1, p2])
+ ptemp = _cyl2cart(ptemp)
+ p1cart = ptemp[0]
+ p2cart = ptemp[1]
+ dpcart = p2cart - p1cart
+
+ # set up components
+ rleft = left_edge[:,0]
+ rright = right_edge[:,0]
+ zleft = left_edge[:,1]
+ zright = right_edge[:,1]
+
+ a = (dpcart[0]**2) + (dpcart[1]**2)
+ b = (2*dpcart[0]*p1cart[0]) + (2*dpcart[1]*p1cart[1])
+ cleft = ((p1cart[0]**2) + (p1cart[1]**2)) - rleft**2
+ cright = ((p1cart[0]**2) + (p1cart[1]**2)) - rright**2
+ twoa = 2*a
+ b2 = b**2
+
+ # Compute positive and negative times and associated masks
+ I = left_edge.shape[0]
+ tmleft = np.empty(I, dtype='float64')
+ tpleft = np.empty(I, dtype='float64')
+ tmright = np.empty(I, dtype='float64')
+ tpright = np.empty(I, dtype='float64')
+ for i in range(I):
+ tmleft[i] = (-b - math.sqrt(b2 - 4*a*cleft[i])) / twoa
+ tpleft[i] = (-b + math.sqrt(b2 - 4*a*cleft[i])) / twoa
+ tmright[i] = (-b - math.sqrt(b2 - 4*a*cright[i])) / twoa
+ tpright[i] = (-b + math.sqrt(b2 - 4*a*cright[i])) / twoa
+
+ tmmright = np.logical_and(~np.isnan(tmright), rright <= p1[0])
+ tpmright = np.logical_and(~np.isnan(tpright), rright <= p2[0])
+
+ tmmleft = np.logical_and(~np.isnan(tmleft), rleft <= p1[0])
+ tpmleft = np.logical_and(~np.isnan(tpleft), rleft <= p2[0])
+
+ # compute first cut of indexes and thetas, which
+ # have been filtered by those values for which intersection
+ # times are impossible (see above masks). Note that this is
+ # still independnent of z.
+ inds = np.unique(np.concatenate([np.argwhere(tmmleft).flat,
+ np.argwhere(tpmleft).flat,
+ np.argwhere(tmmright).flat,
+ np.argwhere(tpmright).flat,]))
+ if 0 == inds.shape[0]:
+ inds = np.arange(I)
+ thetaleft = np.empty(I)
+ thetaleft.fill(p1[2])
+ thetaright = np.empty(I)
+ thetaleft.fill(p2[2])
+ else:
+ thetaleft = np.arctan2((p1cart[1] + tmleft[inds]*dpcart[1]),
+ (p1cart[0] + tmleft[inds]*dpcart[0]))
+ nans = np.isnan(thetaleft)
+ thetaleft[nans] = np.arctan2((p1cart[1] + tpleft[inds[nans]]*dpcart[1]),
+ (p1cart[0] + tpleft[inds[nans]]*dpcart[0]))
+
+ thetaright = np.arctan2((p1cart[1] + tmright[inds]*dpcart[1]),
+ (p1cart[0] + tmright[inds]*dpcart[0]))
+ nans = np.isnan(thetaright)
+ thetaright[nans] = np.arctan2((p1cart[1] + tpright[inds[nans]]*dpcart[1]),
+ (p1cart[0] + tpright[inds[nans]]*dpcart[0]))
+
+
"""\
-D = F - E
-
-Ecart = na.array((E[0]*na.cos(E[2]), E[0]*na.sin(E[2]), E[1]))
-Fcart = na.array((F[0]*na.cos(F[2]), F[0]*na.sin(F[2]), F[1]))
-Dcart = Fcart - Ecart
-
-# <codecell>
-
-rleft = pf.h.grid_left_edge[:,0]
-rright = pf.h.grid_right_edge[:,0]
-zleft = pf.h.grid_left_edge[:,1]
-zright = pf.h.grid_right_edge[:,1]
-
-a = (Dcart**2)[:2].sum()
-b = (2*Dcart*Ecart)[:2].sum()
-cleft = (Ecart**2)[:2].sum() - rleft**2
-cright = (Ecart**2)[:2].sum() - rright**2
-
-tmleft = (-b - na.sqrt(b**2 - 4*a*cleft)) / (2*a)
-tpleft = (-b + na.sqrt(b**2 - 4*a*cleft)) / (2*a)
-tmright = (-b - na.sqrt(b**2 - 4*a*cright)) / (2*a)
-tpright = (-b + na.sqrt(b**2 - 4*a*cright)) / (2*a)
-
-tmmright = np.logical_and(~np.isnan(tmright), rright <= E[0])
-tpmright = np.logical_and(~np.isnan(tpright), rright <= F[0])
-
-tmmleft = np.logical_and(~np.isnan(tmleft), rleft <= E[0])
-tpmleft = np.logical_and(~np.isnan(tpleft), rleft <= F[0])
-
-# <codecell>
-
-ind = np.unique(np.concatenate([np.argwhere(tmmleft).flat, np.argwhere(tpmleft).flat, np.argwhere(tmmright).flat, np.argwhere(tpmright).flat,]))
-
-thetaleft = np.arctan2((Ecart[1] + tmleft[ind]*Dcart[1]), (Ecart[0] + tmleft[ind]*Dcart[0]))
-nans = np.isnan(thetaleft)
-thetaleft[nans] = np.arctan2((Ecart[1] + tpleft[ind[nans]]*Dcart[1]), (Ecart[0] + tpleft[ind[nans]]*Dcart[0]))
-
-thetaright = np.arctan2((Ecart[1] + tmright[ind]*Dcart[1]), (Ecart[0] + tmright[ind]*Dcart[0]))
-nans = np.isnan(thetaleft)
-thetaright[nans] = np.arctan2((Ecart[1] + tpright[ind[nans]]*Dcart[1]), (Ecart[0] + tpright[ind[nans]]*Dcart[0]))
-
-#thetaleft += np.pi*3/2
-#thetaright += np.pi*3/2
-
-if 0 == len(ind):
- print "Ind len zero"
- I = len(rleft)
- ind = np.arange(I)
- thetaleft = np.empty(I)
- thetaleft.fill(E[2])
- thetaright = np.empty(I)
- thetaleft.fill(F[2])
# <codecell>
https://bitbucket.org/yt_analysis/yt-3.0/changeset/864d721e3d1a/
changeset: 864d721e3d1a
branch: yt-3.0
user: scopatz
date: 2012-09-05 22:14:17
summary: Stopping work to make movie for Don.
affected #: 1 file
diff -r 2020eb46759bd8aa14da171111edad63c6aa1469 -r 864d721e3d1a4efd8c27970263f261513a319a9f yt/utilities/lib/alt_ray_tracers.pyx
--- a/yt/utilities/lib/alt_ray_tracers.pyx
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -112,7 +112,7 @@
cleft, cright, thetaleft, thetaright,
tmleft, tpleft, tmright, tpright
cdef np.ndarray[np.int32_t, ndim=1] inds
- cdef np.ndarray[np.float64_t, ndim=2] rztheta, ptemp
+ cdef np.ndarray[np.float64_t, ndim=2] rztheta, ptemp, b1, b2
# set up points
dp = p2 - p1
@@ -180,21 +180,8 @@
thetaright[nans] = np.arctan2((p1cart[1] + tpright[inds[nans]]*dpcart[1]),
(p1cart[0] + tpright[inds[nans]]*dpcart[0]))
-
-"""\
-
-# <codecell>
-
-a = E
-b = F
-
-#c = np.array([rleft[ind], zleft[ind], thetaleft]).T
-#d = np.array([rleft[ind], zright[ind], thetaleft]).T
-
-#c = np.array([rright[ind], zleft[ind], thetaright]).T
-#d = np.array([rright[ind], zright[ind], thetaright]).T
-
-c = np.array([rleft[ind], zright[ind], thetaleft]).T
+ # Set up the cell boundary arrays
+ b1 = np.concatenate() np.array([rleft[ind], zright[ind], thetaleft]).T
c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
c = np.append(c, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
@@ -221,6 +208,7 @@
ind = np.append(ind, origind)
ind = np.append(ind, origind)
+"""\
tsec, intsec = intersect(a, b, c, d)
print tsec
print np.isnan(tsec).all(), len(tsec)
https://bitbucket.org/yt_analysis/yt-3.0/changeset/17ff254e4016/
changeset: 17ff254e4016
branch: yt-3.0
user: scopatz
date: 2012-09-06 04:10:07
summary: cylindrical ray tracer is working.
affected #: 1 file
diff -r 864d721e3d1a4efd8c27970263f261513a319a9f -r 17ff254e4016de55cfeceec099e45212b86e88c4 yt/utilities/lib/alt_ray_tracers.pyx
--- a/yt/utilities/lib/alt_ray_tracers.pyx
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -77,8 +77,7 @@
def clyindrical_ray_trace(np.ndarray[np.float64_t, ndim=1] p1,
np.ndarray[np.float64_t, ndim=1] p2,
np.ndarray[np.float64_t, ndim=2] left_edges,
- np.ndarray[np.float64_t, ndim=2] right_edges)
-
+ np.ndarray[np.float64_t, ndim=2] right_edges):
"""Computes straight (cartesian) rays being traced through a
cylindrical geometry.
@@ -105,14 +104,14 @@
indexes into the grid cells which the ray crosses in order.
"""
- int i, I
- cdef np.float64_t a, b,b2, twoa
- cdef np.ndarray[np.float64_t, ndim=1] dp, p1cart, p2cart, dpcart, t, s,
- rleft, rright, zleft, zright,
- cleft, cright, thetaleft, thetaright,
- tmleft, tpleft, tmright, tpright
- cdef np.ndarray[np.int32_t, ndim=1] inds
- cdef np.ndarray[np.float64_t, ndim=2] rztheta, ptemp, b1, b2
+ cdef int i, I
+ cdef np.float64_t a, b, bsqrd, twoa
+ cdef np.ndarray[np.float64_t, ndim=1] dp, p1cart, p2cart, dpcart, t, s, \
+ rleft, rright, zleft, zright, \
+ cleft, cright, thetaleft, thetaright, \
+ tmleft, tpleft, tmright, tpright, tsect
+ cdef np.ndarray[np.int64_t, ndim=1] inds, tinds, sinds
+ cdef np.ndarray[np.float64_t, ndim=2] xyz, rztheta, ptemp, b1, b2, dsect
# set up points
dp = p2 - p1
@@ -123,29 +122,29 @@
dpcart = p2cart - p1cart
# set up components
- rleft = left_edge[:,0]
- rright = right_edge[:,0]
- zleft = left_edge[:,1]
- zright = right_edge[:,1]
+ rleft = left_edges[:,0]
+ rright = right_edges[:,0]
+ zleft = left_edges[:,1]
+ zright = right_edges[:,1]
a = (dpcart[0]**2) + (dpcart[1]**2)
b = (2*dpcart[0]*p1cart[0]) + (2*dpcart[1]*p1cart[1])
cleft = ((p1cart[0]**2) + (p1cart[1]**2)) - rleft**2
cright = ((p1cart[0]**2) + (p1cart[1]**2)) - rright**2
twoa = 2*a
- b2 = b**2
+ bsqrd = b**2
# Compute positive and negative times and associated masks
- I = left_edge.shape[0]
+ I = left_edges.shape[0]
tmleft = np.empty(I, dtype='float64')
tpleft = np.empty(I, dtype='float64')
tmright = np.empty(I, dtype='float64')
tpright = np.empty(I, dtype='float64')
for i in range(I):
- tmleft[i] = (-b - math.sqrt(b2 - 4*a*cleft[i])) / twoa
- tpleft[i] = (-b + math.sqrt(b2 - 4*a*cleft[i])) / twoa
- tmright[i] = (-b - math.sqrt(b2 - 4*a*cright[i])) / twoa
- tpright[i] = (-b + math.sqrt(b2 - 4*a*cright[i])) / twoa
+ tmleft[i] = (-b - math.sqrt(bsqrd - 4*a*cleft[i])) / twoa
+ tpleft[i] = (-b + math.sqrt(bsqrd - 4*a*cleft[i])) / twoa
+ tmright[i] = (-b - math.sqrt(bsqrd - 4*a*cright[i])) / twoa
+ tpright[i] = (-b + math.sqrt(bsqrd - 4*a*cright[i])) / twoa
tmmright = np.logical_and(~np.isnan(tmright), rright <= p1[0])
tpmright = np.logical_and(~np.isnan(tpright), rright <= p2[0])
@@ -168,12 +167,17 @@
thetaright = np.empty(I)
thetaleft.fill(p2[2])
else:
+ rleft = rleft[inds]
+ rright = rright[inds]
+
+ zleft = zleft[inds]
+ zright = zright[inds]
+
thetaleft = np.arctan2((p1cart[1] + tmleft[inds]*dpcart[1]),
(p1cart[0] + tmleft[inds]*dpcart[0]))
nans = np.isnan(thetaleft)
thetaleft[nans] = np.arctan2((p1cart[1] + tpleft[inds[nans]]*dpcart[1]),
(p1cart[0] + tpleft[inds[nans]]*dpcart[0]))
-
thetaright = np.arctan2((p1cart[1] + tmright[inds]*dpcart[1]),
(p1cart[0] + tmright[inds]*dpcart[0]))
nans = np.isnan(thetaright)
@@ -181,56 +185,46 @@
(p1cart[0] + tpright[inds[nans]]*dpcart[0]))
# Set up the cell boundary arrays
- b1 = np.concatenate() np.array([rleft[ind], zright[ind], thetaleft]).T
-c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
-c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
-c = np.append(c, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
-c = np.append(c, np.array([rleft[ind], zleft[ind], thetaleft]).T, axis=0)
-c = np.append(c, np.array([rright[ind], zright[ind], thetaleft]).T, axis=0)
-c = np.append(c, np.array([rleft[ind], zright[ind], thetaleft]).T, axis=0)
-c = np.append(c, np.array([rright[ind], zleft[ind], thetaleft]).T, axis=0)
+ b1 = np.concatenate([[rleft, zright, thetaleft],
+ [rleft, zleft, thetaleft],
+ [rleft, zleft, thetaleft],
+ [rright, zleft, thetaright],
+ [rleft, zleft, thetaleft],
+ [rright, zright, thetaleft],
+ [rleft, zright, thetaleft],
+ [rright, zleft, thetaleft],
+ ], axis=1).T
-d = np.array([rright[ind], zright[ind], thetaright]).T
-d = np.append(d, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
-d = np.append(d, np.array([rleft[ind], zright[ind], thetaleft]).T, axis=0)
-d = np.append(d, np.array([rright[ind], zright[ind], thetaright]).T, axis=0)
-d = np.append(d, np.array([rleft[ind], zleft[ind], thetaright]).T, axis=0)
-d = np.append(d, np.array([rright[ind], zright[ind], thetaright]).T, axis=0)
-d = np.append(d, np.array([rleft[ind], zright[ind], thetaright]).T, axis=0)
-d = np.append(d, np.array([rright[ind], zleft[ind], thetaright]).T, axis=0)
+ b2 = np.concatenate([[rright, zright, thetaright],
+ [rright, zleft, thetaright],
+ [rleft, zright, thetaleft],
+ [rright, zright, thetaright],
+ [rleft, zleft, thetaright],
+ [rright, zright, thetaright],
+ [rleft, zright, thetaright],
+ [rright, zleft, thetaright],
+ ], axis=1).T
-origind = ind
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
-ind = np.append(ind, origind)
+ inds = np.concatenate([inds, inds, inds, inds, inds, inds, inds, inds])
-"""\
-tsec, intsec = intersect(a, b, c, d)
-print tsec
-print np.isnan(tsec).all(), len(tsec)
-tmask = np.logical_and(0.0<=tsec, tsec<=1.0)
-tsec, utmask = np.unique(tsec[tmask], return_index=True)
-print tsec.max(), tsec.ptp(), len(utmask)
-ind = ind[tmask][utmask]
-xyz = intsec[tmask][utmask]
-s = np.sqrt(((xyz - Ecart)**2).sum(axis=1))
-s, smask = np.unique(s, return_index=True)
-ind = ind[smask]
-xyz = xyz[smask]
-t = s/np.sqrt((Dcart**2).sum())
-si = s.argsort()
-s = s[si]
-t = t[si]
-ind = ind[si]
-xyz = xyz[si]
-rztheta = np.array([np.sqrt(xyz[:,0]**2 + xyz[:,1]**2), xyz[:,2], np.arctan2(xyz[:,1], xyz[:,0])]).T
-rztheta[:,2] = 0.0 + (rztheta[:,2] - np.pi*3/2)%(2*np.pi)
-
-
-
-
-"""
+ # find intersections and compute return values
+ tsect, dsect = _cart_intersect(p1cart, p2cart, _cyl2cart(b1), _cyl2cart(b2))
+ tmask = np.logical_and(0.0<=tsect, tsect<=1.0)
+ tsect, tinds = np.unique(tsect[tmask], return_index=True)
+ inds = inds[tmask][tinds]
+ xyz = dsect[tmask][tinds]
+ s = np.sqrt(((xyz - p1cart)**2).sum(axis=1))
+ s, sinds = np.unique(s, return_index=True)
+ inds = inds[sinds]
+ xyz = xyz[sinds]
+ t = s/np.sqrt((dpcart**2).sum())
+ sinds = s.argsort()
+ s = s[sinds]
+ t = t[sinds]
+ inds = inds[sinds]
+ xyz = xyz[sinds]
+ rztheta = np.concatenate([np.sqrt(xyz[:,0]**2 + xyz[:,1]**2)[:,np.newaxis],
+ xyz[:,2:3],
+ np.arctan2(xyz[:,1], xyz[:,0])[:,np.newaxis]], axis=1)
+ return t, s, rztheta, inds
+ #rztheta[:,2] = 0.0 + (rztheta[:,2] - np.pi*3/2)%(2*np.pi)
https://bitbucket.org/yt_analysis/yt-3.0/changeset/17a8b96a6b39/
changeset: 17a8b96a6b39
branch: yt-3.0
user: scopatz
date: 2012-09-12 00:52:53
summary: tried to integrate cylindrical ray tracing, but is failing for some reason.
affected #: 5 files
diff -r 17ff254e4016de55cfeceec099e45212b86e88c4 -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 yt/data_objects/api.py
--- a/yt/data_objects/api.py
+++ b/yt/data_objects/api.py
@@ -1,4 +1,4 @@
-"""
+"""
API for yt.data_objects
Author: Matthew Turk <matthewturk at gmail.com>
diff -r 17ff254e4016de55cfeceec099e45212b86e88c4 -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -1,4 +1,4 @@
-"""
+"""
The base classes for selecting and returning data.
Author: Matthew Turk <matthewturk at gmail.com>
@@ -61,7 +61,7 @@
return np.zeros(shape, dtype='bool')
def restore_field_information_state(func):
- """
+ """
A decorator that takes a function with the API of (self, grid, field)
and ensures that after the function is called, the field_parameters will
be returned to normal.
@@ -75,14 +75,14 @@
return save_state
class YTFieldData(dict):
- """
+ """
A Container object for field data, instead of just having it be a dict.
"""
pass
class YTDataContainer(object):
- """
+ """
Generic YTDataContainer container. By itself, will attempt to
generate field, read fields (method defined by derived classes)
and deal with passing back and forth field parameters.
@@ -100,7 +100,7 @@
data_object_registry[cls._type_name] = cls
def __init__(self, pf, field_parameters):
- """
+ """
Typically this is never called directly, but only due to inheritance.
It associates a :class:`~yt.data_objects.api.StaticOutput` with the class,
sets its initial set of fields, and the remainder of the arguments
@@ -518,7 +518,7 @@
class YTSelectionContainer2D(YTSelectionContainer):
_key_fields = ['px','py','pdx','pdy']
- """
+ """
Class to represent a set of :class:`YTDataContainer` that's 2-D in nature, and
thus does not have as many actions as the 3-D data types.
"""
diff -r 17ff254e4016de55cfeceec099e45212b86e88c4 -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -1,4 +1,4 @@
-"""
+"""
Data containers based on geometric selection
Author: Matthew Turk <matthewturk at gmail.com>
@@ -35,6 +35,7 @@
from yt.utilities.lib import \
VoxelTraversal, planar_points_in_volume, find_grids_in_inclined_box, \
grid_points_in_volume
+from yt.utilities.lib.alt_ray_tracers import clyindrical_ray_trace
from yt.utilities.orientation import Orientation
from .data_containers import \
YTSelectionContainer1D, YTSelectionContainer2D, YTSelectionContainer3D
@@ -52,8 +53,8 @@
_key_fields = ['x','y','z','dx','dy','dz']
_type_name = "ortho_ray"
_con_args = ('axis', 'coords')
- def __init__(self, axis, coords, pf=None, field_parameters = None):
- """
+ def __init__(self, axis, coords, pf=None, field_parameters=None):
+ """
This is an orthogonal ray cast through the entire domain, at a specific
coordinate.
@@ -103,7 +104,7 @@
_con_args = ('start_point', 'end_point')
sort_by = 't'
def __init__(self, start_point, end_point, pf=None, field_parameters = None):
- """
+ """
This is an arbitrarily-aligned ray cast through the entire domain, at a
specific coordinate.
@@ -141,26 +142,49 @@
#self.vec /= np.sqrt(np.dot(self.vec, self.vec))
self._set_center(self.start_point)
self.set_field_parameter('center', self.start_point)
- self._dts, self._ts = {}, {}
+ self._dts, self._ts, self._masks = {}, {}, {}
def _get_data_from_grid(self, grid, field):
- mask = np.logical_and(self._get_cut_mask(grid),
- grid.child_mask)
- if field == 'dts': return self._dts[grid.id][mask]
- if field == 't': return self._ts[grid.id][mask]
+ if self.pf.geometry == "cylindrical":
+ if grid.id in self._masks:
+ mask = self._masks[grid.id]
+ else:
+ mask = self._get_cut_mask(grid)
+ ts, dts = self._ts[grid.id], self._dts[grid.id]
+ else:
+ mask = np.logical_and(self._get_cut_mask(grid), grid.child_mask)
+ ts, dts = self._ts[grid.id][mask], self._dts[grid.id][mask]
+
+ if field == 'dts':
+ return dts
+ if field == 't':
+ return ts
+
gf = grid[field]
if not iterable(gf):
gf = gf * np.ones(grid.child_mask.shape)
return gf[mask]
def _get_cut_mask(self, grid):
- mask = np.zeros(grid.ActiveDimensions, dtype='int')
- dts = np.zeros(grid.ActiveDimensions, dtype='float64')
- ts = np.zeros(grid.ActiveDimensions, dtype='float64')
- VoxelTraversal(mask, ts, dts, grid.LeftEdge, grid.RightEdge,
- grid.dds, self.center, self.vec)
- self._dts[grid.id] = np.abs(dts)
- self._ts[grid.id] = np.abs(ts)
+ if self.pf.geometry == "cylindrical":
+ _ = clyindrical_ray_trace(self.start_point, self.end_point,
+ grid.LeftEdge, grid.RightEdge)
+ ts, s, rzt, mask = _
+ dts = np.empty(ts.shape, dtype='float64')
+ dts[0], dts[1:] = 0.0, ts[1:] - ts[:-1]
+ grid['r'], grid['z'], grid['theta'] = rzt[:,0], rzt[:,1], rzt[:,2]
+ grid['s'] = s
+ else:
+ mask = np.zeros(grid.ActiveDimensions, dtype='int')
+ dts = np.zeros(grid.ActiveDimensions, dtype='float64')
+ ts = np.zeros(grid.ActiveDimensions, dtype='float64')
+ VoxelTraversal(mask, ts, dts, grid.LeftEdge, grid.RightEdge,
+ grid.dds, self.center, self.vec)
+ dts = np.abs(dts)
+ ts = np.abs(ts)
+ self._dts[grid.id] = dts
+ self._ts[grid.id] = ts
+ self._masks[grid.id] = masks
return mask
diff -r 17ff254e4016de55cfeceec099e45212b86e88c4 -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 yt/utilities/lib/RayIntegrators.pyx
--- a/yt/utilities/lib/RayIntegrators.pyx
+++ b/yt/utilities/lib/RayIntegrators.pyx
@@ -1,4 +1,4 @@
-"""
+"""
Simle integrators for the radiative transfer equation
Author: Matthew Turk <matthewturk at gmail.com>
diff -r 17ff254e4016de55cfeceec099e45212b86e88c4 -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 yt/visualization/plot_collection.py
--- a/yt/visualization/plot_collection.py
+++ b/yt/visualization/plot_collection.py
@@ -1,4 +1,4 @@
-"""
+"""
All of the base-level stuff for plotting.
Author: Matthew Turk <matthewturk at gmail.com>
https://bitbucket.org/yt_analysis/yt-3.0/changeset/81b3e1ff9e39/
changeset: 81b3e1ff9e39
branch: yt-3.0
user: scopatz
date: 2012-09-12 21:47:52
summary: Added RaySelector, not sure if it works.
affected #: 3 files
diff -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 -r 81b3e1ff9e39b106bff7ae8e9ec1798cc931e7f9 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -103,7 +103,7 @@
_type_name = "ray"
_con_args = ('start_point', 'end_point')
sort_by = 't'
- def __init__(self, start_point, end_point, pf=None, field_parameters = None):
+ def __init__(self, start_point, end_point, pf=None, field_parameters=None):
"""
This is an arbitrarily-aligned ray cast through the entire domain, at a
specific coordinate.
@@ -145,6 +145,7 @@
self._dts, self._ts, self._masks = {}, {}, {}
def _get_data_from_grid(self, grid, field):
+ import pdb; pdb.set_trace()
if self.pf.geometry == "cylindrical":
if grid.id in self._masks:
mask = self._masks[grid.id]
diff -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 -r 81b3e1ff9e39b106bff7ae8e9ec1798cc931e7f9 yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -1,4 +1,4 @@
-"""
+"""
Geometry container base class.
Author: Matthew Turk <matthewturk at gmail.com>
diff -r 17a8b96a6b39a740f2da4d1cafb50b8dd67ac8f4 -r 81b3e1ff9e39b106bff7ae8e9ec1798cc931e7f9 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -1,4 +1,4 @@
-"""
+"""
Geometry selection routines.
Author: Matthew Turk <matthewturk at gmail.com>
@@ -724,6 +724,66 @@
ortho_ray_selector = OrthoRaySelector
+cdef class RaySelector(SelectorObject):
+
+ cdef np.float64_t p1[3]
+ cdef np.float64_t p2[3]
+ cdef np.float64_t vec[3]
+
+ def __init__(self, dobj):
+ cdef int i
+ for i in range(3):
+ self.vec[i] = dobj.vec[i]
+ if (0.0 <= self.vec[0]) and (0.0 <= self.vec[1]) and (0.0 <= self.vec[2]):
+ for i in range(3):
+ self.p1[i] = dobj.start_point[i]
+ self.p2[i] = dobj.end_point[i]
+ else:
+ for i in range(3):
+ self.p1[i] = dobj.end_point[i]
+ self.p2[i] = dobj.start_point[i]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int select_grid(self, np.float64_t left_edge[3],
+ np.float64_t right_edge[3]) nogil:
+ if ((self.p1[0] <= left_edge[0]) and (self.p2[0] > right_edge[0]) and
+ (self.p1[1] <= left_edge[1]) and (self.p2[1] > right_edge[1]) and
+ (self.p1[2] <= left_edge[2]) and (self.p2[2] > right_edge[2])):
+ return 1
+ return 0
+
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3],
+ int eterm[3]) nogil:
+ if ((self.p1[0] <= pos[0] - 0.5*dds[0]) and
+ (self.p2[0] > pos[0] + 0.5*dds[0]) and
+ (self.p1[1] <= pos[1] - 0.5*dds[1]) and
+ (self.p2[1] > pos[1] + 0.5*dds[1]) and
+ (self.p1[2] <= pos[2] - 0.5*dds[2]) and
+ (self.p2[2] > pos[2] + 0.5*dds[2])):
+ return 1
+ return 0
+
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void set_bounds(self,
+ np.float64_t left_edge[3], np.float64_t right_edge[3],
+ np.float64_t dds[3], int ind[3][2], int *check):
+ cdef int i
+ for i in range(3):
+ ind[i][0] = <int> ((self.p1[i] - left_edge[i])/dds[i])
+ ind[i][1] = ind[i][0] + 1
+ check[0] = 0
+
+ray_selector = RaySelector
+
cdef class GridCollectionSelector(SelectorObject):
def __init__(self, dobj):
https://bitbucket.org/yt_analysis/yt-3.0/changeset/7e185d2a815b/
changeset: 7e185d2a815b
branch: yt-3.0
user: scopatz
date: 2012-09-12 22:19:17
summary: Removed pdb
affected #: 1 file
diff -r 81b3e1ff9e39b106bff7ae8e9ec1798cc931e7f9 -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -145,7 +145,6 @@
self._dts, self._ts, self._masks = {}, {}, {}
def _get_data_from_grid(self, grid, field):
- import pdb; pdb.set_trace()
if self.pf.geometry == "cylindrical":
if grid.id in self._masks:
mask = self._masks[grid.id]
https://bitbucket.org/yt_analysis/yt-3.0/changeset/3496f043ab4b/
changeset: 3496f043ab4b
branch: yt-3.0
user: scopatz
date: 2012-09-25 01:35:32
summary: Added amrspace() in testing
affected #: 7 files
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -157,9 +157,9 @@
fpu.skip(f)
if hvals['nboundary'] > 0:
fpu.skip(f, 2)
- self.ngridbound = fpu.read_vector(f, 'i')
+ self.ngridbound = fpu.read_vector(f, 'i').astype("int64")
else:
- self.ngridbound = 0
+ self.ngridbound = np.zeros(hvals['nlevelmax'], dtype='int64')
free_mem = fpu.read_attrs(f, (('free_mem', 5, 'i'), ) )
ordering = fpu.read_vector(f, 'c')
fpu.skip(f, 4)
@@ -182,7 +182,7 @@
f.write(fb.read())
f.seek(0)
mylog.debug("Reading domain AMR % 4i (%0.3e, %0.3e)",
- self.domain_id, self.local_oct_count, self.ngridbound)
+ self.domain_id, self.local_oct_count, self.ngridbound.sum())
def _ng(c, l):
if c < self.amr_header['ncpu']:
ng = self.amr_header['numbl'][l, c]
@@ -192,6 +192,7 @@
return ng
min_level = self.pf.min_level
total = 0
+ nx, ny, nz = (((i-1.0)/2.0) for i in self.amr_header['nx'])
for level in range(self.amr_header['nlevelmax']):
# Easier if do this 1-indexed
for cpu in range(self.amr_header['nboundary'] + self.amr_header['ncpu']):
@@ -201,9 +202,9 @@
ind = fpu.read_vector(f, "I").astype("int64")
fpu.skip(f, 2)
pos = np.empty((ng, 3), dtype='float64')
- pos[:,0] = fpu.read_vector(f, "d")
- pos[:,1] = fpu.read_vector(f, "d")
- pos[:,2] = fpu.read_vector(f, "d")
+ pos[:,0] = fpu.read_vector(f, "d") - nx
+ pos[:,1] = fpu.read_vector(f, "d") - ny
+ pos[:,2] = fpu.read_vector(f, "d") - nz
#pos *= self.pf.domain_width
#pos += self.parameter_file.domain_left_edge
fpu.skip(f, 31)
@@ -327,7 +328,8 @@
def _initialize_oct_handler(self):
self.domains = [RAMSESDomainFile(self.parameter_file, i + 1)
for i in range(self.parameter_file['ncpu'])]
- total_octs = sum(dom.local_oct_count for dom in self.domains)
+ total_octs = sum(dom.local_oct_count #+ dom.ngridbound.sum()
+ for dom in self.domains)
self.num_grids = total_octs
#this merely allocates space for the oct tree
#and nothing else
@@ -337,7 +339,7 @@
self.parameter_file.domain_right_edge)
mylog.debug("Allocating %s octs", total_octs)
self.oct_handler.allocate_domains(
- [dom.local_oct_count + dom.ngridbound
+ [dom.local_oct_count #+ dom.ngridbound.sum()
for dom in self.domains])
#this actually reads every oct and loads it into the octree
for dom in self.domains:
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -53,6 +53,9 @@
cdef np.float64_t DLE[3], DRE[3]
cdef public int nocts
cdef public int max_domain
+ cdef Oct* get(self, ppos)
+ cdef void neighbors(self, Oct *, Oct **)
+ cdef void oct_bounds(self, Oct *, np.float64_t *, np.float64_t *)
cdef class RAMSESOctreeContainer(OctreeContainer):
cdef OctAllocationContainer **domains
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -127,6 +127,39 @@
yield (this.ind, this.local_ind, this.domain)
cur = cur.next
+ cdef void oct_bounds(self, Oct *o, np.float64_t *corner, np.float64_t *size):
+ cdef np.float64_t base_dx[3]
+ cdef int i
+ for i in range(3):
+ size[i] = (self.DRE[i] - self.DLE[i]) / (self.nn[i] << o.level)
+ corner[i] = o.pos[i] * size[i] + self.DLE[i]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef Oct *get(self, ppos):
+ cdef np.int64_t ind[3]
+ cdef np.float64_t dds[3], cp[3], pp[3]
+ cdef Oct *cur
+ cdef int i
+ for i in range(3):
+ pp[i] = ppos[i] - self.DLE[i]
+ dds[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
+ ind[i] = <np.int64_t> (pp[i]/dds[i])
+ cp[i] = (ind[i] + 0.5) * dds[i]
+ cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
+ while cur.children[0][0][0] != NULL:
+ for i in range(3):
+ dds[i] = dds[i] / 2.0
+ if cp[i] > pp[i]:
+ ind[i] = 0
+ cp[i] -= dds[i] / 2.0
+ else:
+ ind[i] = 1
+ cp[i] += dds[i]/2.0
+ cur = cur.children[ind[0]][ind[1]][ind[2]]
+ return cur
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -154,6 +187,93 @@
count[o.domain - 1] += mask[oi,i]
return count
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void neighbors(self, Oct* o, Oct* neighbors[27]):
+ cdef np.int64_t curopos[3]
+ cdef np.int64_t curnpos[3]
+ cdef np.int64_t npos[3]
+ cdef int i, j, k, ni, nj, nk, ind[3], nn, dl, skip
+ cdef np.float64_t dds[3], cp[3], pp[3]
+ cdef Oct* candidate
+ for i in range(27): neighbors[i] = NULL
+ nn = 0
+ for ni in range(3):
+ for nj in range(3):
+ for nk in range(3):
+ #print "Handling neighbor", nn
+ if ni == nj == nk == 1:
+ neighbors[nn] = o
+ nn += 1
+ continue
+ npos[0] = o.pos[0] + (ni - 1)
+ npos[1] = o.pos[1] + (nj - 1)
+ npos[2] = o.pos[2] + (nk - 1)
+ for i in range(3):
+ # Periodicity
+ if npos[i] == -1:
+ npos[i] = (self.nn[i] << o.level) - 1
+ elif npos[i] == (self.nn[i] << o.level):
+ npos[i] = 0
+ curopos[i] = o.pos[i]
+ curnpos[i] = npos[i]
+ # Now we have our neighbor position and a safe place to
+ # keep it. curnpos will be the root index of the neighbor
+ # at a given level, and npos will be constant. curopos is
+ # the candidate root at a level.
+ candidate = o
+ while candidate != NULL:
+ if ((curopos[0] == curnpos[0]) and
+ (curopos[1] == curnpos[1]) and
+ (curopos[2] == curnpos[2])):
+ break
+ # This one doesn't meet it, so we pop up a level.
+ # First we update our positions, then we update our
+ # candidate.
+ for i in range(3):
+ # We strip a digit off the right
+ curopos[i] = (curopos[i] >> 1)
+ curnpos[i] = (curnpos[i] >> 1)
+ # Now we update to the candidate's parent, which should
+ # have a matching position to curopos[]
+ candidate = candidate.parent
+ if candidate == NULL:
+ # Worst case scenario
+ for i in range(3):
+ ind[i] = (npos[i] >> (o.level))
+ #print "NULL", npos[0], npos[1], npos[2], ind[0], ind[1], ind[2], o.level
+ candidate = self.root_mesh[ind[0]][ind[1]][ind[2]]
+ # Now we have the common root, which may be NULL
+ while candidate.level < o.level:
+ dl = o.level - (candidate.level + 1)
+ for i in range(3):
+ ind[i] = (npos[i] >> dl) & 1
+ if candidate.children[0][0][0] == NULL: break
+ candidate = candidate.children[ind[0]][ind[1]][ind[2]]
+ neighbors[nn] = candidate
+ nn += 1
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def get_neighbor_boundaries(self, ppos):
+ cdef Oct *main = self.get(ppos)
+ cdef Oct* neighbors[27]
+ self.neighbors(main, neighbors)
+ cdef np.ndarray[np.float64_t, ndim=2] bounds
+ cdef np.float64_t corner[3], size[3]
+ bounds = np.zeros((27,6), dtype="float64")
+ cdef int i, ii
+ cdef np.float64_t base_dx[3]
+ tnp = 0
+ for i in range(27):
+ self.oct_bounds(neighbors[i], corner, size)
+ for ii in range(3):
+ bounds[i, ii] = corner[ii]
+ bounds[i, 3+ii] = size[ii]
+ return bounds
+
cdef class RAMSESOctreeContainer(OctreeContainer):
def allocate_domains(self, domain_counts):
@@ -215,7 +335,7 @@
@cython.cdivision(True)
def add(self, int curdom, int curlevel, int ng,
np.ndarray[np.float64_t, ndim=2] pos,
- int local_domain):
+ int local_domain, int skip_boundary = 1):
cdef int level, no, p, i, j, k, ind[3]
cdef int local = (local_domain == curdom)
cdef Oct *cur, *next = NULL
@@ -224,15 +344,20 @@
if curdom > self.max_domain: curdom = local_domain
cdef OctAllocationContainer *cont = self.domains[curdom - 1]
cdef int initial = cont.n_assigned
+ cdef int in_boundary = 0
# How do we bootstrap ourselves?
for p in range(no):
#for every oct we're trying to add find the
#floating point unitary position on this level
+ in_boundary = 0
for i in range(3):
pp[i] = pos[p, i]
dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
ind[i] = <np.int64_t> (pp[i]/dds[i])
cp[i] = (ind[i] + 0.5) * dds[i]
+ if ind[i] < 0 or ind[i] >= self.nn[i]:
+ in_boundary = 1
+ if skip_boundary == in_boundary == 1: continue
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
if curlevel != 0:
@@ -422,6 +547,20 @@
cdef ParticleArrays *last_sd
cdef Oct** oct_list
+ def allocate_root(self):
+ cdef int i, j, k
+ cdef Oct *cur
+ for i in range(self.nn[0]):
+ for j in range(self.nn[1]):
+ for k in range(self.nn[2]):
+ cur = self.allocate_oct()
+ cur.level = 0
+ cur.pos[0] = i
+ cur.pos[1] = j
+ cur.pos[2] = k
+ cur.parent = NULL
+ self.root_mesh[i][j][k] = cur
+
def __dealloc__(self):
cdef i, j, k
for i in range(self.nn[0]):
@@ -444,8 +583,89 @@
free(o.sd.domain_id)
free(o)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def icoords(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ np.int64_t cell_count):
+ cdef np.ndarray[np.int64_t, ndim=2] coords
+ coords = np.empty((cell_count, 3), dtype="int64")
+ cdef int oi, i, ci, ii
+ ci = 0
+ for oi in range(self.nocts):
+ o = self.oct_list[oi]
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ ii = ((k*2)+j)*2+i
+ if mask[oi, ii] == 1:
+ coords[ci, 0] = (o.pos[0] << 1) + i
+ coords[ci, 1] = (o.pos[1] << 1) + j
+ coords[ci, 2] = (o.pos[2] << 1) + k
+ ci += 1
+ return coords
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def ires(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ np.int64_t cell_count):
+ cdef np.ndarray[np.int64_t, ndim=1] res
+ res = np.empty(cell_count, dtype="int64")
+ cdef int oi, i, ci
+ ci = 0
+ for oi in range(self.nocts):
+ o = self.oct_list[oi]
+ for i in range(8):
+ if mask[oi, i] == 1:
+ res[ci] = o.level
+ ci += 1
+ return res
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def fcoords(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ np.int64_t cell_count):
+ cdef np.ndarray[np.float64_t, ndim=2] coords
+ coords = np.empty((cell_count, 3), dtype="float64")
+ cdef int oi, i, ci
+ cdef np.float64_t base_dx[3], dx[3], pos[3]
+ for i in range(3):
+ # This is the base_dx, but not the base distance from the center
+ # position. Note that the positions will also all be offset by
+ # dx/2.0. This is also for *oct grids*, not cells.
+ base_dx[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
+ ci = 0
+ cdef int proc
+ for oi in range(self.nocts):
+ proc = 0
+ for i in range(8):
+ if mask[oi, i] == 1:
+ proc = 1
+ break
+ if proc == 0: continue
+ o = self.oct_list[oi]
+ for i in range(3):
+ # This gives the *grid* width for this level
+ dx[i] = base_dx[i] / (1 << o.level)
+ # o.pos is the *grid* index, so pos[i] is the center of the
+ # first cell in the grid
+ pos[i] = self.DLE[i] + o.pos[i]*dx[i] + dx[i]/4.0
+ dx[i] = dx[i] / 2.0 # This is now the *offset*
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ ii = ((k*2)+j)*2+i
+ if mask[oi, ii] == 0: continue
+ coords[ci, 0] = pos[0] + dx[0] * i
+ coords[ci, 1] = pos[1] + dx[1] * j
+ coords[ci, 2] = pos[2] + dx[2] * k
+ ci += 1
+ return coords
+
def allocate_domains(self, domain_counts):
- cdef int count, i
+ pass
def finalize(self):
self.oct_list = <Oct**> malloc(sizeof(Oct*)*self.nocts)
@@ -509,6 +729,7 @@
cdef np.float64_t dds[3], cp[3], pp[3]
cdef int ind[3]
self.max_domain = max(self.max_domain, domain_id)
+ if self.root_mesh[0][0][0] == NULL: self.allocate_root()
for p in range(no):
level = 0
for i in range(3):
@@ -518,10 +739,7 @@
cp[i] = (ind[i] + 0.5) * dds[i]
cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
if cur == NULL:
- cur = self.allocate_oct()
- self.root_mesh[ind[0]][ind[1]][ind[2]] = cur
- for i in range(3):
- cur.pos[i] = ind[i] # root level
+ raise RuntimeError
if cur.sd.np == 32:
self.refine_oct(cur, cp)
while cur.sd.np < 0:
@@ -539,6 +757,7 @@
self.refine_oct(cur, cp)
# Now we copy in our particle
pi = cur.sd.np
+ cur.level = level
for i in range(3):
cur.sd.pos[i][pi] = pp[i]
cur.sd.domain_id[pi] = domain_id
@@ -551,9 +770,11 @@
for j in range(2):
for k in range(2):
noct = self.allocate_oct()
- noct.pos[0] = o.pos[0] << 1 + i
- noct.pos[1] = o.pos[1] << 1 + j
- noct.pos[2] = o.pos[2] << 1 + k
+ noct.level = o.level + 1
+ noct.pos[0] = (o.pos[0] << 1) + i
+ noct.pos[1] = (o.pos[1] << 1) + j
+ noct.pos[2] = (o.pos[2] << 1) + k
+ noct.parent = o
o.children[i][j][k] = noct
for m in range(32):
for i in range(3):
@@ -615,3 +836,18 @@
for i in range(o.sd.np):
dmask[o.sd.domain_id[i]] = 1
return dmask.astype("bool")
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def count_neighbor_particles(self, ppos):
+ cdef Oct *main = self.get(ppos)
+ cdef Oct* neighbors[27]
+ self.neighbors(main, neighbors)
+ cdef int i, ni, dl, tnp
+ tnp = 0
+ for i in range(27):
+ if neighbors[i].sd != NULL:
+ tnp += neighbors[i].sd.np
+ return tnp
+
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/geometry/setup.py
--- a/yt/geometry/setup.py
+++ b/yt/geometry/setup.py
@@ -9,13 +9,19 @@
config = Configuration('geometry',parent_package,top_path)
config.add_extension("oct_container",
["yt/geometry/oct_container.pyx"],
- libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ libraries=["m"],
+ depends=["yt/utilities/lib/fp_utils.pxd",
+ "yt/geometry/oct_container.pxd",
+ "yt/geometry/selection_routines.pxd"])
config.add_extension("selection_routines",
["yt/geometry/selection_routines.pyx"],
extra_compile_args=['-fopenmp'],
extra_link_args=['-fopenmp'],
include_dirs=["yt/utilities/lib/"],
- libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ libraries=["m"],
+ depends=["yt/utilities/lib/fp_utils.pxd",
+ "yt/geometry/oct_container.pxd",
+ "yt/geometry/selection_routines.pxd"])
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
return config
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/testing.py
--- /dev/null
+++ b/yt/testing.py
@@ -0,0 +1,127 @@
+"""Provides utility and helper functions for testing in yt.
+
+Author: Anthony Scpatz <scopatz at gmail.com>
+Affiliation: The University of Chicago
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Anthony Scopatz. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+
+def amrspace(extent, levels=7, cells=8):
+ """Creates two numpy arrays representing the left and right bounds of
+ an AMR grid as well as an array for the AMR level of each cell.
+
+ Parameters
+ ----------
+ extent : array-like
+ This a sequence of length 2*ndims that is the bounds of each dimension.
+ For example, the 2D unit square would be given by [0.0, 1.0, 0.0, 1.0].
+ A 3D cylindrical grid may look like [0.0, 2.0, -1.0, 1.0, 0.0, 2*np.pi].
+ levels : int or sequence of ints, optional
+ This is the number of AMR refinement levels. If given as a sequence (of
+ length ndims), then each dimension will be refined down to this level.
+ All values in this array must be the same or zero. A zero valued dimension
+ indicates that this dim should not be refined. Taking the 3D cylindrical
+ example above if we don't want refine theta but want r and z at 5 we would
+ set levels=(5, 5, 0).
+ cells : int, optional
+ This is the number of cells per refinement level.
+
+ Returns
+ -------
+ left : float ndarray, shape=(npoints, ndims)
+ The left AMR grid points.
+ right : float ndarray, shape=(npoints, ndims)
+ The right AMR grid points.
+ level : int ndarray, shape=(npoints,)
+ The AMR level for each point.
+
+ Examples
+ --------
+ >>> l, r, lvl = amrspace([0.0, 2.0, 1.0, 2.0, 0.0, 3.14], levels=(3,3,0), cells=2)
+ >>> print l
+ [[ 0. 1. 0. ]
+ [ 0.25 1. 0. ]
+ [ 0. 1.125 0. ]
+ [ 0.25 1.125 0. ]
+ [ 0.5 1. 0. ]
+ [ 0. 1.25 0. ]
+ [ 0.5 1.25 0. ]
+ [ 1. 1. 0. ]
+ [ 0. 1.5 0. ]
+ [ 1. 1.5 0. ]]
+
+ """
+ extent = np.asarray(extent, dtype='f8')
+ dextent = extent[1::2] - extent[::2]
+ ndims = len(dextent)
+
+ if isinstance(levels, int):
+ minlvl = maxlvl = levels
+ levels = np.array([levels]*ndims, dtype='int32')
+ else:
+ levels = np.asarray(levels, dtype='int32')
+ minlvl = levels.min()
+ maxlvl = levels.max()
+ if minlvl != maxlvl and (minlvl != 0 or set([minlvl, maxlvl]) != set(levels)):
+ raise ValueError("all levels must have the same value or zero.")
+ dims_zero = (levels == 0)
+ dims_nonzero = ~dims_zero
+ ndims_nonzero = dims_nonzero.sum()
+
+ npoints = (cells**ndims_nonzero - 1)*maxlvl + 1
+ left = np.empty((npoints, ndims), dtype='float64')
+ right = np.empty((npoints, ndims), dtype='float64')
+ level = np.empty(npoints, dtype='int32')
+
+ # fill zero dims
+ left[:,dims_zero] = extent[::2][dims_zero]
+ right[:,dims_zero] = extent[1::2][dims_zero]
+
+ # fill non-zero dims
+ dcell = 1.0 / cells
+ left_slice = tuple([slice(extent[2*n], extent[2*n+1], extent[2*n+1]) if \
+ dims_zero[n] else slice(0.0,1.0,dcell) for n in range(ndims)])
+ right_slice = tuple([slice(extent[2*n+1], extent[2*n], -extent[2*n+1]) if \
+ dims_zero[n] else slice(dcell,1.0+dcell,dcell) for n in range(ndims)])
+ left_norm_grid = np.reshape(np.mgrid[left_slice].T.flat[ndims:], (-1, ndims))
+ lng_zero = left_norm_grid[:,dims_zero]
+ lng_nonzero = left_norm_grid[:,dims_nonzero]
+
+ right_norm_grid = np.reshape(np.mgrid[right_slice].T.flat[ndims:], (-1, ndims))
+ rng_zero = right_norm_grid[:,dims_zero]
+ rng_nonzero = right_norm_grid[:,dims_nonzero]
+
+ level[0] = maxlvl
+ left[0,:] = extent[::2]
+ right[0,dims_zero] = extent[1::2][dims_zero]
+ right[0,dims_nonzero] = (dcell**maxlvl)*dextent[dims_nonzero] + extent[::2][dims_nonzero]
+ for i, lvl in enumerate(range(maxlvl, 0, -1)):
+ start = (cells**ndims_nonzero - 1)*i + 1
+ stop = (cells**ndims_nonzero - 1)*(i+1) + 1
+ dsize = dcell**(lvl-1) * dextent[dims_nonzero]
+ level[start:stop] = lvl
+ left[start:stop,dims_zero] = lng_zero
+ left[start:stop,dims_nonzero] = lng_nonzero*dsize + extent[::2][dims_nonzero]
+ right[start:stop,dims_zero] = rng_zero
+ right[start:stop,dims_nonzero] = rng_nonzero*dsize + extent[::2][dims_nonzero]
+
+ return left, right, level
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/utilities/linear_interpolators.py
--- a/yt/utilities/linear_interpolators.py
+++ b/yt/utilities/linear_interpolators.py
@@ -53,7 +53,8 @@
my_vals = np.zeros(x_vals.shape, dtype='float64')
lib.UnilinearlyInterpolate(self.table, x_vals, self.x_bins, x_i, my_vals)
- return my_vals.reshape(orig_shape)
+ my_vals.shape = orig_shape
+ return my_vals
class BilinearFieldInterpolator:
def __init__(self, table, boundaries, field_names, truncate=False):
@@ -86,7 +87,8 @@
lib.BilinearlyInterpolate(self.table,
x_vals, y_vals, self.x_bins, self.y_bins,
x_i, y_i, my_vals)
- return my_vals.reshape(orig_shape)
+ my_vals.shape = orig_shape
+ return my_vals
class TrilinearFieldInterpolator:
def __init__(self, table, boundaries, field_names, truncate = False):
@@ -125,31 +127,8 @@
x_vals, y_vals, z_vals,
self.x_bins, self.y_bins, self.z_bins,
x_i, y_i, z_i, my_vals)
- return my_vals.reshape(orig_shape)
-
- # Use notation from Paul Bourke's page on interpolation
- # http://local.wasp.uwa.edu.au/~pbourke/other/interpolation/
- x = (x_vals - self.x_bins[x_i]) / (self.x_bins[x_i+1] - self.x_bins[x_i])
- y = (y_vals - self.y_bins[y_i]) / (self.y_bins[y_i+1] - self.y_bins[y_i])
- z = (z_vals - self.z_bins[z_i]) / (self.z_bins[z_i+1] - self.z_bins[z_i])
- xm = (self.x_bins[x_i+1] - x_vals) / (self.x_bins[x_i+1] - self.x_bins[x_i])
- ym = (self.y_bins[y_i+1] - y_vals) / (self.y_bins[y_i+1] - self.y_bins[y_i])
- zm = (self.z_bins[z_i+1] - z_vals) / (self.z_bins[z_i+1] - self.z_bins[z_i])
- if np.any(np.isnan(self.table)):
- raise ValueError
- if np.any(np.isnan(x) | np.isnan(y) | np.isnan(z)):
- raise ValueError
- if np.any(np.isnan(xm) | np.isnan(ym) | np.isnan(zm)):
- raise ValueError
- my_vals = self.table[x_i ,y_i ,z_i ] * (xm*ym*zm)
- my_vals += self.table[x_i+1,y_i ,z_i ] * (x *ym*zm)
- my_vals += self.table[x_i ,y_i+1,z_i ] * (xm*y *zm)
- my_vals += self.table[x_i ,y_i ,z_i+1] * (xm*ym*z )
- my_vals += self.table[x_i+1,y_i ,z_i+1] * (x *ym*z )
- my_vals += self.table[x_i ,y_i+1,z_i+1] * (xm*y *z )
- my_vals += self.table[x_i+1,y_i+1,z_i ] * (x *y *zm)
- my_vals += self.table[x_i+1,y_i+1,z_i+1] * (x *y *z )
- return my_vals.reshape(orig_shape)
+ my_vals.shape = orig_shape
+ return my_vals
def get_centers(pf, filename, center_cols, radius_col, unit='1'):
"""
diff -r 7e185d2a815b79f07b2ea164b5763393a74c0c8e -r 3496f043ab4b3b809e73963c150c7afc561674f1 yt/utilities/tests/test_interpolators.py
--- /dev/null
+++ b/yt/utilities/tests/test_interpolators.py
@@ -0,0 +1,13 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+import yt.utilities.linear_interpolators as lin
+
+def setup():
+ pass
+
+
+def test_linear_interpolator():
+ random_data = np.random.random(128)
+ x = {"Random":np.mgrid[0.0:1.0:128j]}
+ ufi = lin.UnilinearFieldInterpolator(random_data, (0.0, 1.0), "Random", True)
+ assert_array_equal(ufi(x), random_data)
https://bitbucket.org/yt_analysis/yt-3.0/changeset/7844f5a2ee9e/
changeset: 7844f5a2ee9e
branch: yt-3.0
user: scopatz
date: 2012-09-25 03:59:32
summary: First cut at ray trace tests.
affected #: 1 file
diff -r 3496f043ab4b3b809e73963c150c7afc561674f1 -r 7844f5a2ee9e04627e7234418419e40be401ed9e yt/utilities/lib/tests/test_alt_ray_tracers.py
--- /dev/null
+++ b/yt/utilities/lib/tests/test_alt_ray_tracers.py
@@ -0,0 +1,60 @@
+"""Tests for non-cartesian ray tracers."""
+import nose
+import numpy as np
+
+from nose.tools import assert_equal, assert_not_equal, assert_raises, raises, \
+ assert_almost_equal, assert_true, assert_false, assert_in
+from numpy.testing import assert_array_equal, assert_array_almost_equal
+from yt.testing import amrspace
+
+from yt.utilities.lib.alt_ray_tracers import clyindrical_ray_trace
+
+left_grid = right_grid = amr_levels = center_grid = data = None
+
+def setup():
+ # set up some sample cylindrical grid data, radiating out from center
+ global left_grid, right_grid, amr_levels, center_grid
+ np.seterr(all='ignore')
+ l1, r1, lvl1 = amrspace([0.0, 1.0, 0.0, -1.0, 0.0, 2*np.pi], levels=(7,7,0))
+ l2, r2, lvl2 = amrspace([0.0, 1.0, 0.0, 1.0, 0.0, 2*np.pi], levels=(7,7,0))
+ left_grid = np.concatenate([l1,l2], axis=0)
+ right_grid = np.concatenate([r1,r2], axis=0)
+ amr_levels = np.concatenate([lvl1,lvl2], axis=0)
+ center_grid = (left_grid + right_grid) / 2.0
+ data = np.cos(np.sqrt(np.sum(center_grid[:,:2]**2, axis=1)))**2 # cos^2
+
+
+point_pairs = [
+ # p1 p2
+ ([0.5, -1.0, 0.0], [1.0, 1.0, 0.75*np.pi]), # Everything different
+ ([0.5, -1.0, 0.0], [0.5, 1.0, 0.75*np.pi]), # r same
+ ([0.5, -1.0, 0.0], [0.5, 1.0, np.pi]), # diagonal through z-axis
+ # straight through z-axis
+ ([0.5, 0.0, 0.0], [0.5, 0.0, np.pi]),
+ #([0.5, 0.0, np.pi*3/2 + 0.0], [0.5, 0.0, np.pi*3/2 + np.pi]),
+ #([0.5, 0.0, np.pi/2 + 0.0], [0.5, 0.0, np.pi/2 + np.pi]),
+ #([0.5, 0.0, np.pi + 0.0], [0.5, 0.0, np.pi + np.pi]),
+ # const z, not through z-axis
+ ([0.5, 0.1, 0.0], [0.5, 0.1, 0.75*np.pi]),
+ #([0.5, 0.1, np.pi + 0.0], [0.5, 0.1, np.pi + 0.75*np.pi]),
+ #([0.5, 0.1, np.pi*3/2 + 0.0], [0.5, 0.1, np.pi*3/2 + 0.75*np.pi]),
+ #([0.5, 0.1, np.pi/2 + 0.0], [0.5, 0.1, np.pi/2 + 0.75*np.pi]),
+ #([0.5, 0.1, 2*np.pi + 0.0], [0.5, 0.1, 2*np.pi + 0.75*np.pi]),
+ #([0.5, 0.1, np.pi/4 + 0.0], [0.5, 0.1, np.pi/4 + 0.75*np.pi]),
+ #([0.5, 0.1, np.pi*3/8 + 0.0], [0.5, 0.1, np.pi*3/8 + 0.75*np.pi]),
+ ([0.5, -1.0, 0.75*np.pi], [1.0, 1.0, 0.75*np.pi]), # r,z different - theta same
+ ([0.5, -1.0, 0.75*np.pi], [0.5, 1.0, 0.75*np.pi]), # z-axis parallel
+ ([0.0, -1.0, 0.0], [0.0, 1.0, 0.0]), # z-axis itself
+ ]
+
+
+def check_monotonic_inc(arr):
+ assert_true(np.all(0.0 <= (arr[1:] - arr[:-1])))
+
+
+def test_clyindrical_ray_trace():
+ for p1, p2 in point_pairs:
+ p1, p2 = np.asarray(p1), np.asarray(p2)
+ t, s, rztheta, inds = clyindrical_ray_trace(p1, p2, left_grid, right_grid)
+ yield check_monotonic_inc, t
+ yield check_monotonic_inc, s
https://bitbucket.org/yt_analysis/yt-3.0/changeset/52a1b2ecea94/
changeset: 52a1b2ecea94
branch: yt-3.0
user: scopatz
date: 2012-09-25 04:42:49
summary: finished writing up some tests for cylindrical coord ray tracer.
affected #: 1 file
diff -r 7844f5a2ee9e04627e7234418419e40be401ed9e -r 52a1b2ecea9488edb1df3f0090fcd18bac38586e yt/utilities/lib/tests/test_alt_ray_tracers.py
--- a/yt/utilities/lib/tests/test_alt_ray_tracers.py
+++ b/yt/utilities/lib/tests/test_alt_ray_tracers.py
@@ -3,11 +3,12 @@
import numpy as np
from nose.tools import assert_equal, assert_not_equal, assert_raises, raises, \
- assert_almost_equal, assert_true, assert_false, assert_in
+ assert_almost_equal, assert_true, assert_false, assert_in, assert_less_equal, \
+ assert_greater_equal
from numpy.testing import assert_array_equal, assert_array_almost_equal
from yt.testing import amrspace
-from yt.utilities.lib.alt_ray_tracers import clyindrical_ray_trace
+from yt.utilities.lib.alt_ray_tracers import clyindrical_ray_trace, _cyl2cart
left_grid = right_grid = amr_levels = center_grid = data = None
@@ -24,7 +25,7 @@
data = np.cos(np.sqrt(np.sum(center_grid[:,:2]**2, axis=1)))**2 # cos^2
-point_pairs = [
+point_pairs = np.array([
# p1 p2
([0.5, -1.0, 0.0], [1.0, 1.0, 0.75*np.pi]), # Everything different
([0.5, -1.0, 0.0], [0.5, 1.0, 0.75*np.pi]), # r same
@@ -45,16 +46,40 @@
([0.5, -1.0, 0.75*np.pi], [1.0, 1.0, 0.75*np.pi]), # r,z different - theta same
([0.5, -1.0, 0.75*np.pi], [0.5, 1.0, 0.75*np.pi]), # z-axis parallel
([0.0, -1.0, 0.0], [0.0, 1.0, 0.0]), # z-axis itself
- ]
+ ])
def check_monotonic_inc(arr):
assert_true(np.all(0.0 <= (arr[1:] - arr[:-1])))
+def check_bounds(arr, blower, bupper):
+ assert_true(np.all(blower <= arr))
+ assert_true(np.all(bupper >= arr))
+
def test_clyindrical_ray_trace():
- for p1, p2 in point_pairs:
- p1, p2 = np.asarray(p1), np.asarray(p2)
+ for pair in point_pairs:
+ p1, p2 = pair
+ p1cart, p2cart = _cyl2cart(pair)
+ pathlen = np.sqrt(np.sum((p2cart - p1cart)**2))
+
t, s, rztheta, inds = clyindrical_ray_trace(p1, p2, left_grid, right_grid)
+ npoints = len(t)
+
yield check_monotonic_inc, t
+ yield assert_less_equal, 0.0, t[0]
+ yield assert_less_equal, t[-1], 1.0
+
yield check_monotonic_inc, s
+ yield assert_less_equal, 0.0, s[0]
+ yield assert_less_equal, s[-1], pathlen
+ yield assert_equal, npoints, len(s)
+
+ yield assert_equal, (npoints, 3), rztheta.shape
+ yield check_bounds, rztheta[:,0], 0.0, 1.0
+ yield check_bounds, rztheta[:,1], -1.0, 1.0
+ yield check_bounds, rztheta[:,2], 0.0, 2*np.pi
+ yield check_monotonic_inc, rztheta[:,2]
+
+ yield assert_equal, npoints, len(inds)
+ yield check_bounds, inds, 0, len(left_grid)-1
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
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