[Yt-svn] commit/yt: 8 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Oct 6 12:57:36 PDT 2011
8 new changesets in yt:
http://bitbucket.org/yt_analysis/yt/changeset/4b1015184456/
changeset: 4b1015184456
branch: yt
user: MatthewTurk
date: 2011-10-03 21:50:17
summary: Adding data_coords option to the image line callback, and adding support for
phase plots to same.
affected #: 1 file (-1 bytes)
--- a/yt/visualization/plot_modifications.py Mon Oct 03 10:34:29 2011 -0700
+++ b/yt/visualization/plot_modifications.py Mon Oct 03 15:50:17 2011 -0400
@@ -50,8 +50,14 @@
pass
def convert_to_pixels(self, plot, coord, offset = True):
- x0, x1 = plot.xlim
- y0, y1 = plot.ylim
+ if plot.xlim is not None:
+ x0, x1 = plot.xlim
+ else:
+ x0, x1 = plot._axes.get_xlim()
+ if plot.ylim is not None:
+ y0, y1 = plot.ylim
+ else:
+ y0, y1 = plot._axes.get_ylim()
l, b, width, height = mpl_get_bounds(plot._axes.bbox)
dx = width / (x1-x0)
dy = height / (y1-y0)
@@ -473,7 +479,7 @@
class ImageLineCallback(LinePlotCallback):
_type_name = "image_line"
- def __init__(self, p1, p2, plot_args = None):
+ def __init__(self, p1, p2, data_coords=False, plot_args = None):
"""
Plot from *p1* to *p2* (image plane coordinates)
with *plot_args* fed into the plot.
@@ -484,19 +490,27 @@
if plot_args is None: plot_args = {}
self.plot_args = plot_args
self._ids = []
+ self.data_coords = data_coords
def __call__(self, plot):
# We manually clear out any previous calls to this callback:
plot._axes.lines = [l for l in plot._axes.lines if id(l) not in self._ids]
- p1 = self.convert_to_pixels(plot, self.p1)
- p2 = self.convert_to_pixels(plot, self.p2)
+ kwargs = self.plot_args.copy()
+ if self.data_coords and len(plot.image._A.shape) == 2:
+ p1 = self.convert_to_pixels(plot, self.p1)
+ p2 = self.convert_to_pixels(plot, self.p2)
+ else:
+ p1, p2 = self.p1, self.p2
+ if not self.data_coords:
+ kwargs["transform"] = plot._axes.transAxes
+
px, py = (p1[0], p2[0]), (p1[1], p2[1])
# Save state
xx0, xx1 = plot._axes.get_xlim()
yy0, yy1 = plot._axes.get_ylim()
plot._axes.hold(True)
- ii = plot._axes.plot(px, py, **self.plot_args)
+ ii = plot._axes.plot(px, py, **kwargs)
self._ids.append(id(ii[0]))
# Reset state
plot._axes.set_xlim(xx0,xx1)
@@ -905,7 +919,7 @@
class TextLabelCallback(PlotCallback):
_type_name = "text"
- def __init__(self, pos, text, data_coords=False,text_args = None):
+ def __init__(self, pos, text, data_coords=False, text_args = None):
"""
Accepts a position in (0..1, 0..1) of the image, some text and
optionally some text arguments. If data_coords is True,
http://bitbucket.org/yt_analysis/yt/changeset/1c5a9cf00283/
changeset: 1c5a9cf00283
branch: yt
user: MatthewTurk
date: 2011-10-04 16:59:08
summary: Refactoring marching cubes code to be external to partitioned grid.
affected #: 1 file (-1 bytes)
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Mon Oct 03 15:50:17 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 10:59:08 2011 -0400
@@ -70,6 +70,11 @@
Triangle *next
np.float64_t p[3][3]
+cdef struct TriangleCollection:
+ int count
+ Triangle *first
+ Triangle *current
+
cdef Triangle *AddTriangle(Triangle *self,
np.float64_t p0[3], np.float64_t p1[3], np.float64_t p2[3]):
cdef Triangle *nn = <Triangle *> malloc(sizeof(Triangle))
@@ -843,315 +848,25 @@
for i in range(3):
vel[i] /= vel_mag[0]
- #@cython.boundscheck(False)
- #@cython.wraparound(False)
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
@cython.cdivision(True)
def get_isocontour_triangles(self, np.float64_t isovalue, int field_id = 0):
# Much of this was inspired by code from Paul Bourke's website:
# http://paulbourke.net/geometry/polygonise/
- cdef int *edge_table=[
- 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
- 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
- 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
- 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
- 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
- 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
- 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
- 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
- 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
- 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
- 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
- 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
- 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
- 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
- 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
- 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
- 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
- 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
- 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
- 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
- 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
- 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
- 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
- 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
- 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
- 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
- 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
- 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
- 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
- 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
- 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
- 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 ]
-
- cdef int **tri_table = \
- [[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
- [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
- [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
- [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
- [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
- [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
- [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
- [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
- [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
- [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
- [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
- [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
- [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
- [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
- [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
- [5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1],
- [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1],
- [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
- [0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
- [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1],
- [10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1],
- [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1],
- [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1],
- [5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1],
- [9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1],
- [0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1],
- [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1],
- [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1],
- [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1],
- [2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1],
- [7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1],
- [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1],
- [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1],
- [11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1],
- [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1],
- [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1],
- [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1],
- [11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
- [1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1],
- [9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1],
- [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1],
- [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
- [0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
- [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1],
- [6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1],
- [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1],
- [6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1],
- [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1],
- [1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
- [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1],
- [6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1],
- [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1],
- [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
- [3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
- [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1],
- [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1],
- [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1],
- [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1],
- [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1],
- [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
- [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1],
- [10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1],
- [10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1],
- [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1],
- [1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1],
- [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1],
- [0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1],
- [10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1],
- [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1],
- [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1],
- [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1],
- [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
- [3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1],
- [6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1],
- [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1],
- [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1],
- [10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1],
- [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
- [7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1],
- [7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1],
- [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1],
- [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1],
- [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1],
- [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
- [0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1],
- [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
- [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
- [2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
- [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1],
- [7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1],
- [2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1],
- [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1],
- [10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1],
- [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1],
- [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1],
- [7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1],
- [6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1],
- [8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1],
- [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1],
- [6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1],
- [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1],
- [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1],
- [8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1],
- [0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1],
- [1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1],
- [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1],
- [10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1],
- [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
- [10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
- [5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
- [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1],
- [9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
- [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1],
- [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1],
- [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
- [7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1],
- [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1],
- [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1],
- [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
- [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1],
- [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
- [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1],
- [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1],
- [6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1],
- [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1],
- [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1],
- [6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1],
- [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1],
- [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1],
- [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1],
- [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1],
- [9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1],
- [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
- [1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
- [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1],
- [0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1],
- [5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1],
- [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1],
- [11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1],
- [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1],
- [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
- [2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1],
- [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1],
- [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1],
- [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
- [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1],
- [9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1],
- [9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1],
- [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1],
- [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1],
- [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1],
- [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1],
- [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1],
- [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
- [9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1],
- [5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1],
- [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
- [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1],
- [8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1],
- [0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1],
- [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1],
- [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1],
- [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1],
- [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1],
- [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1],
- [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1],
- [11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1],
- [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1],
- [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1],
- [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
- [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1],
- [1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1],
- [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1],
- [4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1],
- [0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1],
- [3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1],
- [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1],
- [0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1],
- [9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1],
- [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
- [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]]
+ # Cython makes us toss this in here, which I think will change in a
+ # future release.
cdef int i, j, k, n
cdef int offset
cdef np.float64_t gv[8]
cdef int cubeindex
- cdef np.float64_t vertlist[12][3]
- cdef int ntriang = 0
cdef np.float64_t *intdata = NULL
cdef np.float64_t x, y, z
cdef np.float64_t mu
- cdef Triangle *first = NULL
- cdef Triangle *current = NULL
+ cdef TriangleCollection triangles
+ triangles.first = triangles.current = NULL
+ triangles.count = 0
x = self.left_edge[0]
for i in range(self.dims[0]):
y = self.left_edge[1]
@@ -1162,68 +877,373 @@
+ j * (self.dims[2] + 1) + k
intdata = self.data[field_id] + offset
offset_fill(self.dims, intdata, gv)
- cubeindex = 0
- for n in range(8):
- if gv[n] < isovalue:
- cubeindex |= (1 << n)
- if edge_table[cubeindex] == 0:
- z += self.dds[2]
- continue
- if (edge_table[cubeindex] & 1): # 0,0,0 with 1,0,0
- vertex_interp(gv[0], gv[1], isovalue, vertlist[0],
- self.dds, x, y, z, 0, 1)
- if (edge_table[cubeindex] & 2): # 1,0,0 with 1,1,0
- vertex_interp(gv[1], gv[2], isovalue, vertlist[1],
- self.dds, x, y, z, 1, 2)
- if (edge_table[cubeindex] & 4): # 1,1,0 with 0,1,0
- vertex_interp(gv[2], gv[3], isovalue, vertlist[2],
- self.dds, x, y, z, 2, 3)
- if (edge_table[cubeindex] & 8): # 0,1,0 with 0,0,0
- vertex_interp(gv[3], gv[0], isovalue, vertlist[3],
- self.dds, x, y, z, 3, 0)
- if (edge_table[cubeindex] & 16): # 0,0,1 with 1,0,1
- vertex_interp(gv[4], gv[5], isovalue, vertlist[4],
- self.dds, x, y, z, 4, 5)
- if (edge_table[cubeindex] & 32): # 1,0,1 with 1,1,1
- vertex_interp(gv[5], gv[6], isovalue, vertlist[5],
- self.dds, x, y, z, 5, 6)
- if (edge_table[cubeindex] & 64): # 1,1,1 with 0,1,1
- vertex_interp(gv[6], gv[7], isovalue, vertlist[6],
- self.dds, x, y, z, 6, 7)
- if (edge_table[cubeindex] & 128): # 0,1,1 with 0,0,1
- vertex_interp(gv[7], gv[4], isovalue, vertlist[7],
- self.dds, x, y, z, 7, 4)
- if (edge_table[cubeindex] & 256): # 0,0,0 with 0,0,1
- vertex_interp(gv[0], gv[4], isovalue, vertlist[8],
- self.dds, x, y, z, 0, 4)
- if (edge_table[cubeindex] & 512): # 1,0,0 with 1,0,1
- vertex_interp(gv[1], gv[5], isovalue, vertlist[9],
- self.dds, x, y, z, 1, 5)
- if (edge_table[cubeindex] & 1024): # 1,1,0 with 1,1,1
- vertex_interp(gv[2], gv[6], isovalue, vertlist[10],
- self.dds, x, y, z, 2, 6)
- if (edge_table[cubeindex] & 2048): # 0,1,0 with 0,1,1
- vertex_interp(gv[3], gv[7], isovalue, vertlist[11],
- self.dds, x, y, z, 3, 7)
- n = 0
- while 1:
- current = AddTriangle(current,
- vertlist[tri_table[cubeindex][n ]],
- vertlist[tri_table[cubeindex][n+1]],
- vertlist[tri_table[cubeindex][n+2]])
- ntriang += 1
- if first == NULL: first = current
- n += 3
- if tri_table[cubeindex][n] == -1: break
+ march_cubes(gv, isovalue, self.dds, x, y, z,
+ &triangles)
z += self.dds[2]
y += self.dds[1]
x += self.dds[0]
# Hallo, we are all done.
cdef np.ndarray[np.float64_t, ndim=2] vertices
- vertices = np.zeros((ntriang*3,3), dtype='float64')
- FillAndWipeTriangles(vertices, first)
+ vertices = np.zeros((triangles.count*3,3), dtype='float64')
+ print "Found ", triangles.count, "triangles"
+ FillAndWipeTriangles(vertices, triangles.first)
return vertices
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef march_cubes(np.float64_t gv[8], np.float64_t isovalue,
+ np.float64_t dds[3],
+ np.float64_t x, np.float64_t y, np.float64_t z,
+ TriangleCollection *triangles):
+ cdef int *edge_table=[
+ 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
+ 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+ 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
+ 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+ 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
+ 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+ 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
+ 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+ 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
+ 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+ 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
+ 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+ 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
+ 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+ 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
+ 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+ 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
+ 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+ 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
+ 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+ 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
+ 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+ 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
+ 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
+ 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
+ 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
+ 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
+ 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
+ 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
+ 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
+ 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
+ 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 ]
+
+ cdef int **tri_table = \
+ [[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
+ [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
+ [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
+ [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
+ [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
+ [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
+ [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
+ [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
+ [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
+ [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
+ [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
+ [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
+ [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
+ [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
+ [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
+ [5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1],
+ [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1],
+ [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
+ [0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
+ [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1],
+ [10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1],
+ [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1],
+ [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1],
+ [5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1],
+ [9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1],
+ [0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1],
+ [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1],
+ [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1],
+ [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1],
+ [2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1],
+ [7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1],
+ [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1],
+ [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1],
+ [11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1],
+ [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1],
+ [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1],
+ [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1],
+ [11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
+ [1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1],
+ [9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1],
+ [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1],
+ [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
+ [0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
+ [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1],
+ [6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1],
+ [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1],
+ [6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1],
+ [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1],
+ [1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
+ [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1],
+ [6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1],
+ [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1],
+ [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
+ [3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
+ [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1],
+ [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1],
+ [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1],
+ [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1],
+ [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1],
+ [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
+ [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1],
+ [10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1],
+ [10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1],
+ [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1],
+ [1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1],
+ [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1],
+ [0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1],
+ [10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1],
+ [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1],
+ [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1],
+ [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1],
+ [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
+ [3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1],
+ [6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1],
+ [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1],
+ [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1],
+ [10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1],
+ [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
+ [7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1],
+ [7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1],
+ [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1],
+ [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1],
+ [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1],
+ [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
+ [0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1],
+ [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
+ [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
+ [2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
+ [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1],
+ [7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1],
+ [2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1],
+ [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1],
+ [10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1],
+ [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1],
+ [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1],
+ [7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1],
+ [6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1],
+ [8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1],
+ [6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1],
+ [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1],
+ [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1],
+ [8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1],
+ [0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1],
+ [1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1],
+ [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1],
+ [10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1],
+ [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
+ [10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
+ [5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
+ [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1],
+ [9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
+ [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1],
+ [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1],
+ [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
+ [7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1],
+ [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1],
+ [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1],
+ [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
+ [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1],
+ [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
+ [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1],
+ [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1],
+ [6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1],
+ [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1],
+ [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1],
+ [6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1],
+ [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1],
+ [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1],
+ [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1],
+ [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1],
+ [9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1],
+ [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
+ [1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
+ [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1],
+ [0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1],
+ [5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1],
+ [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1],
+ [11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1],
+ [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1],
+ [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
+ [2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1],
+ [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1],
+ [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1],
+ [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
+ [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1],
+ [9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1],
+ [9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1],
+ [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1],
+ [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1],
+ [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1],
+ [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1],
+ [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1],
+ [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
+ [9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1],
+ [5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1],
+ [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
+ [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1],
+ [8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1],
+ [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1],
+ [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1],
+ [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1],
+ [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1],
+ [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1],
+ [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1],
+ [11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1],
+ [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1],
+ [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1],
+ [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
+ [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1],
+ [1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1],
+ [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1],
+ [4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1],
+ [0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1],
+ [3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1],
+ [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1],
+ [0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1],
+ [9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1],
+ [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+ [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]]
+ cdef np.float64_t vertlist[12][3]
+ cdef int cubeindex = 0
+ cdef int n
+ for n in range(8):
+ if gv[n] < isovalue:
+ cubeindex |= (1 << n)
+ if edge_table[cubeindex] == 0: return
+ if (edge_table[cubeindex] & 1): # 0,0,0 with 1,0,0
+ vertex_interp(gv[0], gv[1], isovalue, vertlist[0],
+ dds, x, y, z, 0, 1)
+ if (edge_table[cubeindex] & 2): # 1,0,0 with 1,1,0
+ vertex_interp(gv[1], gv[2], isovalue, vertlist[1],
+ dds, x, y, z, 1, 2)
+ if (edge_table[cubeindex] & 4): # 1,1,0 with 0,1,0
+ vertex_interp(gv[2], gv[3], isovalue, vertlist[2],
+ dds, x, y, z, 2, 3)
+ if (edge_table[cubeindex] & 8): # 0,1,0 with 0,0,0
+ vertex_interp(gv[3], gv[0], isovalue, vertlist[3],
+ dds, x, y, z, 3, 0)
+ if (edge_table[cubeindex] & 16): # 0,0,1 with 1,0,1
+ vertex_interp(gv[4], gv[5], isovalue, vertlist[4],
+ dds, x, y, z, 4, 5)
+ if (edge_table[cubeindex] & 32): # 1,0,1 with 1,1,1
+ vertex_interp(gv[5], gv[6], isovalue, vertlist[5],
+ dds, x, y, z, 5, 6)
+ if (edge_table[cubeindex] & 64): # 1,1,1 with 0,1,1
+ vertex_interp(gv[6], gv[7], isovalue, vertlist[6],
+ dds, x, y, z, 6, 7)
+ if (edge_table[cubeindex] & 128): # 0,1,1 with 0,0,1
+ vertex_interp(gv[7], gv[4], isovalue, vertlist[7],
+ dds, x, y, z, 7, 4)
+ if (edge_table[cubeindex] & 256): # 0,0,0 with 0,0,1
+ vertex_interp(gv[0], gv[4], isovalue, vertlist[8],
+ dds, x, y, z, 0, 4)
+ if (edge_table[cubeindex] & 512): # 1,0,0 with 1,0,1
+ vertex_interp(gv[1], gv[5], isovalue, vertlist[9],
+ dds, x, y, z, 1, 5)
+ if (edge_table[cubeindex] & 1024): # 1,1,0 with 1,1,1
+ vertex_interp(gv[2], gv[6], isovalue, vertlist[10],
+ dds, x, y, z, 2, 6)
+ if (edge_table[cubeindex] & 2048): # 0,1,0 with 0,1,1
+ vertex_interp(gv[3], gv[7], isovalue, vertlist[11],
+ dds, x, y, z, 3, 7)
+ n = 0
+ while 1:
+ print "Adding triangle"
+ triangles.current = AddTriangle(triangles.current,
+ vertlist[tri_table[cubeindex][n ]],
+ vertlist[tri_table[cubeindex][n+1]],
+ vertlist[tri_table[cubeindex][n+2]])
+ triangles.count += 1
+ if triangles.first == NULL:
+ triangles.first = triangles.current
+ n += 3
+ if tri_table[cubeindex][n] == -1: break
+
+
cdef class GridFace:
cdef int direction
cdef public np.float64_t coord
http://bitbucket.org/yt_analysis/yt/changeset/87e2f9f6f53b/
changeset: 87e2f9f6f53b
branch: yt
user: MatthewTurk
date: 2011-10-04 19:22:23
summary: Adding method to walk over grids instead of partitioned bricks to get
isocontours.
affected #: 1 file (-1 bytes)
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 10:59:08 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 13:22:23 2011 -0400
@@ -885,13 +885,54 @@
# Hallo, we are all done.
cdef np.ndarray[np.float64_t, ndim=2] vertices
vertices = np.zeros((triangles.count*3,3), dtype='float64')
- print "Found ", triangles.count, "triangles"
FillAndWipeTriangles(vertices, triangles.first)
return vertices
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+def march_cubes_grid(np.float64_t isovalue,
+ np.ndarray[np.float64_t, ndim=3] values,
+ np.ndarray[np.int32_t, ndim=3] mask,
+ np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] dxs):
+ cdef int dims[3]
+ cdef int i, j, k, n
+ cdef int offset
+ cdef np.float64_t gv[8]
+ cdef np.float64_t *intdata = NULL
+ cdef np.float64_t x, y, z
+ cdef TriangleCollection triangles
+ for i in range(3): dims[i] = values.shape[i] - 1
+ triangles.first = triangles.current = NULL
+ triangles.count = 0
+ cdef np.float64_t *data = <np.float64_t *> values.data
+ cdef np.float64_t *dds = <np.float64_t *> dxs.data
+ x = left_edge[0]
+ for i in range(dims[0]):
+ y = left_edge[1]
+ for j in range(dims[1]):
+ z = left_edge[2]
+ for k in range(dims[2]):
+ if mask[i,j,k] == 1:
+ offset = i * (dims[1] + 1) * (dims[2] + 1) \
+ + j * (dims[2] + 1) + k
+ intdata = data + offset
+ offset_fill(dims, intdata, gv)
+ march_cubes(gv, isovalue, dds, x, y, z,
+ &triangles)
+ z += dds[2]
+ y += dds[1]
+ x += dds[0]
+ # Hallo, we are all done.
+ cdef np.ndarray[np.float64_t, ndim=2] vertices
+ vertices = np.zeros((triangles.count*3,3), dtype='float64')
+ FillAndWipeTriangles(vertices, triangles.first)
+ return vertices
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
cdef march_cubes(np.float64_t gv[8], np.float64_t isovalue,
np.float64_t dds[3],
np.float64_t x, np.float64_t y, np.float64_t z,
@@ -1232,7 +1273,6 @@
dds, x, y, z, 3, 7)
n = 0
while 1:
- print "Adding triangle"
triangles.current = AddTriangle(triangles.current,
vertlist[tri_table[cubeindex][n ]],
vertlist[tri_table[cubeindex][n+1]],
http://bitbucket.org/yt_analysis/yt/changeset/963bd0b00092/
changeset: 963bd0b00092
branch: yt
user: MatthewTurk
date: 2011-10-05 04:00:49
summary: * Adding the ability to march cubes and calculate fluxes.
* Added helper functions to AMR3DData
More to come. Sample code:
http://paste.yt-project.org/show/1846/
affected #: 2 files (-1 bytes)
--- a/yt/data_objects/data_containers.py Tue Oct 04 13:22:23 2011 -0400
+++ b/yt/data_objects/data_containers.py Tue Oct 04 22:00:49 2011 -0400
@@ -40,7 +40,8 @@
from yt.data_objects.particle_io import particle_handler_registry
from yt.utilities.amr_utils import find_grids_in_inclined_box, \
grid_points_in_volume, planar_points_in_volume, VoxelTraversal, \
- QuadTree, get_box_grids_below_level, ghost_zone_interpolate
+ QuadTree, get_box_grids_below_level, ghost_zone_interpolate, \
+ march_cubes_grid, march_cubes_grid_flux
from yt.utilities.data_point_utilities import CombineGrids, \
DataCubeRefine, DataCubeReplace, FillRegion, FillBuffer
from yt.utilities.definitions import axis_names, x_dict, y_dict
@@ -2411,6 +2412,41 @@
__quantities = None
quantities = property(__get_quantities)
+ def extract_isocontours(self, field, value, filename = None,
+ recenter = False):
+ verts = []
+ for g in self._grids:
+ mask = self._get_cut_mask(g) * g.child_mask
+ vals = g.get_vertex_centered_data(field)
+ my_verts = march_cubes_grid(value, vals, mask, g.LeftEdge, g.dds)
+ verts.append(my_verts)
+ verts = na.concatenate(verts)
+ if recenter:
+ mi = na.min(verts, axis=0)
+ ma = na.max(verts, axis=0)
+ verts = (verts - mi) / (ma - mi).max()
+ if filename is not None:
+ f = open(filename, "w")
+ for v1 in verts:
+ f.write("v %0.16e %0.16e %0.16e\n" % (v1[0], v1[1], v1[2]))
+ for i in range(len(verts)/3):
+ f.write("f %s %s %s\n" % (i*3+1, i*3+2, i*3+3))
+ return verts
+
+ def calculate_isocontour_flux(self, field, value,
+ field_x, field_y, field_z):
+ flux = 0.0
+ for g in self._grids:
+ mask = self._get_cut_mask(g) * g.child_mask
+ vals = g.get_vertex_centered_data(field)
+ xv, yv, zv = [g.get_vertex_centered_data(f) for f in
+ [field_x, field_y, field_z]]
+ my_flux = march_cubes_grid_flux(value, vals, xv, yv, zv,
+ mask, g.LeftEdge, g.dds)
+ print my_flux
+ flux += my_flux
+ return flux
+
def extract_connected_sets(self, field, num_levels, min_val, max_val,
log_space=True, cumulative=True, cache=False):
"""
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 13:22:23 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 22:00:49 2011 -0400
@@ -98,6 +98,14 @@
this = this.next
return count
+cdef void WipeTriangles(Triangle *first):
+ cdef Triangle *this = first
+ cdef Triangle *last
+ while this != NULL:
+ last = this
+ this = this.next
+ free(last)
+
cdef void FillAndWipeTriangles(np.ndarray[np.float64_t, ndim=2] vertices,
Triangle *first):
cdef int count = 0
@@ -933,7 +941,126 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef march_cubes(np.float64_t gv[8], np.float64_t isovalue,
+def march_cubes_grid_flux(
+ np.float64_t isovalue,
+ np.ndarray[np.float64_t, ndim=3] values,
+ np.ndarray[np.float64_t, ndim=3] v1,
+ np.ndarray[np.float64_t, ndim=3] v2,
+ np.ndarray[np.float64_t, ndim=3] v3,
+ np.ndarray[np.int32_t, ndim=3] mask,
+ np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] dxs):
+ cdef int dims[3]
+ cdef int i, j, k, n, m
+ cdef int offset
+ cdef np.float64_t gv[8]
+ cdef np.float64_t *intdata = NULL
+ cdef TriangleCollection triangles
+ cdef Triangle *current = NULL
+ cdef Triangle *last = NULL
+ cdef np.float64_t *data = <np.float64_t *> values.data
+ cdef np.float64_t *v1data = <np.float64_t *> v1.data
+ cdef np.float64_t *v2data = <np.float64_t *> v2.data
+ cdef np.float64_t *v3data = <np.float64_t *> v3.data
+ cdef np.float64_t *dds = <np.float64_t *> dxs.data
+ cdef np.float64_t flux = 0.0
+ cdef np.float64_t center[3], point[3], normal[3], wval, temp, area, s
+ cdef np.float64_t cell_pos[3], fv[3], idds[3]
+ for i in range(3):
+ dims[i] = values.shape[i] - 1
+ idds[i] = 1.0 / dds[i]
+ triangles.first = triangles.current = NULL
+ triangles.count = 0
+ cell_pos[0] = left_edge[0]
+ for i in range(dims[0]):
+ cell_pos[1] = left_edge[1]
+ for j in range(dims[1]):
+ cell_pos[2] = left_edge[2]
+ for k in range(dims[2]):
+ if mask[i,j,k] == 1:
+ offset = i * (dims[1] + 1) * (dims[2] + 1) \
+ + j * (dims[2] + 1) + k
+ intdata = data + offset
+ offset_fill(dims, intdata, gv)
+ march_cubes(gv, isovalue, dds,
+ cell_pos[0], cell_pos[1], cell_pos[2],
+ &triangles)
+ # Now our triangles collection has a bunch. We now
+ # calculate fluxes for each.
+ if last == NULL and triangles.first != NULL:
+ current = triangles.first
+ last = NULL
+ elif last != NULL:
+ current = last.next
+ while current != NULL:
+ # Calculate the center of the triangle
+ wval = 0.0
+ for n in range(3):
+ center[n] = 0.0
+ for n in range(3):
+ for m in range(3):
+ point[m] = (current.p[n][m]-cell_pos[m])*idds[m]
+ # Now we calculate the value at this point
+ temp = offset_interpolate(dims, point, intdata)
+ #print "something", temp, point[0], point[1], point[2]
+ wval += temp
+ for m in range(3):
+ center[m] += temp * point[m]
+ # Now we divide by our normalizing factor
+ for n in range(3):
+ center[n] /= wval
+ # We have our center point of the triangle, in 0..1
+ # coordinates. So now we interpolate our three
+ # fields.
+ fv[0] = offset_interpolate(dims, center, v1data + offset)
+ fv[1] = offset_interpolate(dims, center, v2data + offset)
+ fv[2] = offset_interpolate(dims, center, v3data + offset)
+ # We interpolate again the actual value data
+ wval = offset_interpolate(dims, center, intdata)
+ # Now we have our flux vector and our field value!
+ # We just need a cross product with which we can
+ # dot it.
+ normal[0] = ( current.p[0][1] * current.p[1][2]
+ - current.p[0][2] * current.p[1][1])
+ normal[1] = ( current.p[0][2] * current.p[1][0]
+ - current.p[0][0] * current.p[1][2])
+ normal[2] = ( current.p[0][0] * current.p[1][1]
+ - current.p[0][1] * current.p[1][0])
+ # Dump this somewhere for now
+ temp = wval * (fv[0] * normal[0] +
+ fv[1] * normal[1] +
+ fv[2] * normal[2])
+ # Now we need the area of the triangle. This will take
+ # a lot of time to calculate compared to the rest.
+ # We use Heron's formula.
+ for n in range(3):
+ fv[n] = 0.0
+ for n in range(3):
+ fv[0] += (current.p[0][n] - current.p[2][n])**2.0
+ fv[1] += (current.p[1][n] - current.p[0][n])**2.0
+ fv[2] += (current.p[2][n] - current.p[1][n])**2.0
+ s = 0.0
+ for n in range(3):
+ fv[n] = fv[n]**0.5
+ s += 0.5 * fv[n]
+ area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))
+ area = area**0.5
+ flux += temp*area
+ last = current
+ if current.next == NULL: break
+ current = current.next
+ cell_pos[2] += dds[2]
+ cell_pos[1] += dds[1]
+ cell_pos[0] += dds[0]
+ # Hallo, we are all done.
+ WipeTriangles(triangles.first)
+ return flux
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef int march_cubes(
+ np.float64_t gv[8], np.float64_t isovalue,
np.float64_t dds[3],
np.float64_t x, np.float64_t y, np.float64_t z,
TriangleCollection *triangles):
@@ -1231,10 +1358,12 @@
cdef np.float64_t vertlist[12][3]
cdef int cubeindex = 0
cdef int n
+ cdef int nt = 0
for n in range(8):
if gv[n] < isovalue:
cubeindex |= (1 << n)
- if edge_table[cubeindex] == 0: return
+ if edge_table[cubeindex] == 0:
+ return 0
if (edge_table[cubeindex] & 1): # 0,0,0 with 1,0,0
vertex_interp(gv[0], gv[1], isovalue, vertlist[0],
dds, x, y, z, 0, 1)
@@ -1278,10 +1407,12 @@
vertlist[tri_table[cubeindex][n+1]],
vertlist[tri_table[cubeindex][n+2]])
triangles.count += 1
+ nt += 1
if triangles.first == NULL:
triangles.first = triangles.current
n += 3
if tri_table[cubeindex][n] == -1: break
+ return nt
cdef class GridFace:
http://bitbucket.org/yt_analysis/yt/changeset/ba2c2aa3fe55/
changeset: ba2c2aa3fe55
branch: yt
user: MatthewTurk
date: 2011-10-05 04:12:45
summary: Adding docstrings
affected #: 1 file (-1 bytes)
--- a/yt/data_objects/data_containers.py Tue Oct 04 22:00:49 2011 -0400
+++ b/yt/data_objects/data_containers.py Tue Oct 04 22:12:45 2011 -0400
@@ -2413,7 +2413,54 @@
quantities = property(__get_quantities)
def extract_isocontours(self, field, value, filename = None,
- recenter = False):
+ rescale = False):
+ r"""This identifies isocontours on a cell-by-cell basis, with no
+ consideration of global connectedness, and returns the vertices of the
+ Triangles in that isocontour.
+
+ This function simply returns the vertices of all the triangles
+ calculated by the marching cubes algorithm; for more complex
+ operations, such as identifying connected sets of cells above a given
+ threshold, see the extract_connected_sets function. This is more
+ useful for calculating, for instance, total isocontour area, or
+ visualizing in an external program (such as `MeshLab
+ <http://meshlab.sf.net>`_.)
+
+ Parameters
+ ----------
+ field : string
+ Any field that can be obtained in a data object. This is the field
+ which will be isocontoured.
+ value : float
+ The value at which the isocontour should be calculated.
+ filename : string, optional
+ If supplied, this file will be filled with the vertices in .obj
+ format. Suitable for loading into meshlab.
+ rescale : bool, optional
+ If true, the vertices will be rescaled within their min/max.
+
+ Returns
+ -------
+ verts : array of floats
+ The array of vertices, x,y,z. Taken in threes, these are the
+ triangle vertices.
+
+ References
+ ----------
+
+ .. [1] Marching Cubes: http://en.wikipedia.org/wiki/Marching_cubes
+
+ Examples
+ --------
+ This will create a data object, find a nice value in the center, and
+ output the vertices to "triangles.obj" after rescaling them.
+
+ >>> dd = pf.h.all_data()
+ >>> rho = dd.quantities["WeightedAverageQuantity"](
+ ... "Density", weight="CellMassMsun")
+ >>> verts = dd.extract_isocontours("Density", rho,
+ ... "triangles.obj", True)
+ """
verts = []
for g in self._grids:
mask = self._get_cut_mask(g) * g.child_mask
@@ -2421,7 +2468,7 @@
my_verts = march_cubes_grid(value, vals, mask, g.LeftEdge, g.dds)
verts.append(my_verts)
verts = na.concatenate(verts)
- if recenter:
+ if rescale:
mi = na.min(verts, axis=0)
ma = na.max(verts, axis=0)
verts = (verts - mi) / (ma - mi).max()
@@ -2435,16 +2482,65 @@
def calculate_isocontour_flux(self, field, value,
field_x, field_y, field_z):
+ r"""This identifies isocontours on a cell-by-cell basis, with no
+ consideration of global connectedness, and calculates the flux over
+ those contours.
+
+ This function will conduct marching cubes on all the cells in a given
+ data container (grid-by-grid), and then for each identified triangular
+ segment of an isocontour in a given cell, calculate the normal and area
+ of that triangle, the area of that triangle, and then return:
+
+ area * local_value * (n dot v)
+
+ Where area, local_value, and the vector v are interpolated at the barycenter
+ (weighted by the vertex values) of the triangle.
+
+ Parameters
+ ----------
+ field : string
+ Any field that can be obtained in a data object. This is the field
+ which will be isocontoured and used as the "local_value" in the
+ flux equation.
+ value : float
+ The value at which the isocontour should be calculated.
+ field_x : string
+ The x-component field
+ field_y : string
+ The y-component field
+ field_z : string
+ The z-component field
+
+ Returns
+ -------
+ flux : float
+ The summed flux. Note that it is not currently scaled; this is
+ simply the code-unit area times the fields.
+
+ References
+ ----------
+
+ .. [1] Marching Cubes: http://en.wikipedia.org/wiki/Marching_cubes
+
+ Examples
+ --------
+ This will create a data object, find a nice value in the center, and
+ calculate a flux over it.
+
+ >>> dd = pf.h.all_data()
+ >>> rho = dd.quantities["WeightedAverageQuantity"](
+ ... "Density", weight="CellMassMsun")
+ >>> flux = dd.calculate_isocontour_flux("Density", rho,
+ ... "x-velocity", "y-velocity", "z-velocity")
+ """
flux = 0.0
for g in self._grids:
mask = self._get_cut_mask(g) * g.child_mask
vals = g.get_vertex_centered_data(field)
xv, yv, zv = [g.get_vertex_centered_data(f) for f in
[field_x, field_y, field_z]]
- my_flux = march_cubes_grid_flux(value, vals, xv, yv, zv,
+ flux += march_cubes_grid_flux(value, vals, xv, yv, zv,
mask, g.LeftEdge, g.dds)
- print my_flux
- flux += my_flux
return flux
def extract_connected_sets(self, field, num_levels, min_val, max_val,
http://bitbucket.org/yt_analysis/yt/changeset/b88517d95eca/
changeset: b88517d95eca
branch: yt
user: MatthewTurk
date: 2011-10-06 20:28:25
summary: In response to Sam's comments, fix calculation of normal vector. Check test
code at http://paste.yt-project.org/show/1861/ .
affected #: 3 files (-1 bytes)
--- a/yt/utilities/_amr_utils/FixedInterpolator.c Tue Oct 04 22:12:45 2011 -0400
+++ b/yt/utilities/_amr_utils/FixedInterpolator.c Thu Oct 06 14:28:25 2011 -0400
@@ -127,8 +127,8 @@
return vz[0];
}
-npy_float64 eval_gradient(int *ds, int *ci, npy_float64 *dp,
- npy_float64 *data, npy_float64 *grad)
+void eval_gradient(int ds[3], npy_float64 dp[3],
+ npy_float64 *data, npy_float64 grad[3])
{
// We just take some small value
@@ -145,9 +145,9 @@
//fprintf(stderr, "DIM: %d %0.3lf %0.3lf\n", i, plus, minus);
denom = plus - minus;
dp[i] = plus;
- grad[i] += trilinear_interpolate(ds, ci, dp, data) / denom;
+ grad[i] += offset_interpolate(ds, dp, data) / denom;
dp[i] = minus;
- grad[i] -= trilinear_interpolate(ds, ci, dp, data) / denom;
+ grad[i] -= offset_interpolate(ds, dp, data) / denom;
dp[i] = backup;
normval += grad[i]*grad[i];
}
--- a/yt/utilities/_amr_utils/FixedInterpolator.h Tue Oct 04 22:12:45 2011 -0400
+++ b/yt/utilities/_amr_utils/FixedInterpolator.h Thu Oct 06 14:28:25 2011 -0400
@@ -41,8 +41,7 @@
npy_float64 trilinear_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
npy_float64 *data);
-npy_float64 eval_gradient(int ds[3], int ci[3], npy_float64 dp[3],
- npy_float64 *data, npy_float64 *grad);
+void eval_gradient(int ds[3], npy_float64 dp[3], npy_float64 *data, npy_float64 grad[3]);
void vertex_interp(npy_float64 v1, npy_float64 v2, npy_float64 isovalue,
npy_float64 vl[3], npy_float64 dds[3],
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Tue Oct 04 22:12:45 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Oct 06 14:28:25 2011 -0400
@@ -127,8 +127,8 @@
np.float64_t offset_interpolate(int ds[3], np.float64_t dp[3], np.float64_t *data)
np.float64_t trilinear_interpolate(int ds[3], int ci[3], np.float64_t dp[3],
np.float64_t *data)
- np.float64_t eval_gradient(int *ds, int *ci, np.float64_t *dp,
- np.float64_t *data, np.float64_t *grad)
+ void eval_gradient(int ds[3], np.float64_t dp[3], np.float64_t *data,
+ np.float64_t grad[3])
void offset_fill(int *ds, np.float64_t *data, np.float64_t *gridval)
void vertex_interp(np.float64_t v1, np.float64_t v2, np.float64_t isovalue,
np.float64_t vl[3], np.float64_t dds[3],
@@ -964,8 +964,8 @@
cdef np.float64_t *v3data = <np.float64_t *> v3.data
cdef np.float64_t *dds = <np.float64_t *> dxs.data
cdef np.float64_t flux = 0.0
- cdef np.float64_t center[3], point[3], normal[3], wval, temp, area, s
- cdef np.float64_t cell_pos[3], fv[3], idds[3]
+ cdef np.float64_t center[3], point[3], wval, temp, area, s
+ cdef np.float64_t cell_pos[3], fv[3], idds[3], normal[3]
for i in range(3):
dims[i] = values.shape[i] - 1
idds[i] = 1.0 / dds[i]
@@ -1018,18 +1018,18 @@
# We interpolate again the actual value data
wval = offset_interpolate(dims, center, intdata)
# Now we have our flux vector and our field value!
- # We just need a cross product with which we can
- # dot it.
- normal[0] = ( current.p[0][1] * current.p[1][2]
- - current.p[0][2] * current.p[1][1])
- normal[1] = ( current.p[0][2] * current.p[1][0]
- - current.p[0][0] * current.p[1][2])
- normal[2] = ( current.p[0][0] * current.p[1][1]
- - current.p[0][1] * current.p[1][0])
+ # We just need a normal vector with which we can
+ # dot it. The normal should be equal to the gradient
+ # in the center of the triangle, or thereabouts.
+ eval_gradient(dims, center, intdata, normal)
+ temp = 0.0
+ for n in range(3):
+ temp += normal[n]*normal[n]
+ temp = temp**0.5
# Dump this somewhere for now
temp = wval * (fv[0] * normal[0] +
fv[1] * normal[1] +
- fv[2] * normal[2])
+ fv[2] * normal[2])/temp
# Now we need the area of the triangle. This will take
# a lot of time to calculate compared to the rest.
# We use Heron's formula.
http://bitbucket.org/yt_analysis/yt/changeset/c3efa0c49ab8/
changeset: c3efa0c49ab8
branch: yt
user: MatthewTurk
date: 2011-10-06 21:40:38
summary: Reflecting comments from Sam, this changes two things:
* The flux is now defined as (optional) fluxing_field * norm dot vec. If the
fluxing field is not defined, it is assumed to be 1.0.
* The docstrings now reflect this.
In the future, some cleverness might be able to do the 1.0 setting internal to
the Cython code, rather than supplying an additional field.
affected #: 2 files (-1 bytes)
--- a/yt/data_objects/data_containers.py Thu Oct 06 14:28:25 2011 -0400
+++ b/yt/data_objects/data_containers.py Thu Oct 06 15:40:38 2011 -0400
@@ -2481,20 +2481,28 @@
return verts
def calculate_isocontour_flux(self, field, value,
- field_x, field_y, field_z):
+ field_x, field_y, field_z, fluxing_field = None):
r"""This identifies isocontours on a cell-by-cell basis, with no
consideration of global connectedness, and calculates the flux over
those contours.
This function will conduct marching cubes on all the cells in a given
data container (grid-by-grid), and then for each identified triangular
- segment of an isocontour in a given cell, calculate the normal and area
- of that triangle, the area of that triangle, and then return:
+ segment of an isocontour in a given cell, calculate the gradient (i.e.,
+ normal) in the isocontoured field, interpolate the local value of the
+ "fluxing" field, the area of the triangle, and then return:
- area * local_value * (n dot v)
+ area * local_flux_value * (n dot v)
Where area, local_value, and the vector v are interpolated at the barycenter
- (weighted by the vertex values) of the triangle.
+ (weighted by the vertex values) of the triangle. Note that this
+ specifically allows for the field fluxing across the surface to be
+ *different* from the field being contoured. If the fluxing_field is
+ not specified, it is assumed to be 1.0 everywhere, and the raw flux
+ with no local-weighting is returned.
+
+ Additionally, the returned flux is defined as flux *into* the surface,
+ not flux *out of* the surface.
Parameters
----------
@@ -2510,6 +2518,9 @@
The y-component field
field_z : string
The z-component field
+ fluxing_field : string, optional
+ The field whose passage over the surface is of interest. If not
+ specified, assumed to be 1.0 everywhere.
Returns
-------
@@ -2525,22 +2536,26 @@
Examples
--------
This will create a data object, find a nice value in the center, and
- calculate a flux over it.
+ calculate the metal flux over it.
>>> dd = pf.h.all_data()
>>> rho = dd.quantities["WeightedAverageQuantity"](
... "Density", weight="CellMassMsun")
>>> flux = dd.calculate_isocontour_flux("Density", rho,
- ... "x-velocity", "y-velocity", "z-velocity")
+ ... "x-velocity", "y-velocity", "z-velocity", "Metal_Density")
"""
flux = 0.0
for g in self._grids:
mask = self._get_cut_mask(g) * g.child_mask
vals = g.get_vertex_centered_data(field)
+ if fluxing_field is None:
+ ff = na.ones(vals.shape, dtype="float64")
+ else:
+ ff = g.get_vertex_centered_data(fluxing_field)
xv, yv, zv = [g.get_vertex_centered_data(f) for f in
[field_x, field_y, field_z]]
flux += march_cubes_grid_flux(value, vals, xv, yv, zv,
- mask, g.LeftEdge, g.dds)
+ ff, mask, g.LeftEdge, g.dds)
return flux
def extract_connected_sets(self, field, num_levels, min_val, max_val,
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Oct 06 14:28:25 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Oct 06 15:40:38 2011 -0400
@@ -947,6 +947,7 @@
np.ndarray[np.float64_t, ndim=3] v1,
np.ndarray[np.float64_t, ndim=3] v2,
np.ndarray[np.float64_t, ndim=3] v3,
+ np.ndarray[np.float64_t, ndim=3] flux_field,
np.ndarray[np.int32_t, ndim=3] mask,
np.ndarray[np.float64_t, ndim=1] left_edge,
np.ndarray[np.float64_t, ndim=1] dxs):
@@ -962,6 +963,7 @@
cdef np.float64_t *v1data = <np.float64_t *> v1.data
cdef np.float64_t *v2data = <np.float64_t *> v2.data
cdef np.float64_t *v3data = <np.float64_t *> v3.data
+ cdef np.float64_t *fdata = <np.float64_t *> flux_field.data
cdef np.float64_t *dds = <np.float64_t *> dxs.data
cdef np.float64_t flux = 0.0
cdef np.float64_t center[3], point[3], wval, temp, area, s
@@ -1016,7 +1018,7 @@
fv[1] = offset_interpolate(dims, center, v2data + offset)
fv[2] = offset_interpolate(dims, center, v3data + offset)
# We interpolate again the actual value data
- wval = offset_interpolate(dims, center, intdata)
+ wval = offset_interpolate(dims, center, fdata + offset)
# Now we have our flux vector and our field value!
# We just need a normal vector with which we can
# dot it. The normal should be equal to the gradient
http://bitbucket.org/yt_analysis/yt/changeset/7efe36e44c0e/
changeset: 7efe36e44c0e
branch: yt
user: MatthewTurk
date: 2011-10-06 21:54:13
summary: Fixing negative sign.
affected #: 1 file (-1 bytes)
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Oct 06 15:40:38 2011 -0400
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Oct 06 15:54:13 2011 -0400
@@ -1027,7 +1027,8 @@
temp = 0.0
for n in range(3):
temp += normal[n]*normal[n]
- temp = temp**0.5
+ # Take the negative, to ensure it points inwardly
+ temp = -(temp**0.5)
# Dump this somewhere for now
temp = wval * (fv[0] * normal[0] +
fv[1] * normal[1] +
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
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