[yt-svn] commit/yt: 18 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon May 23 20:30:35 PDT 2016
18 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/5d7165b50f9d/
Changeset: 5d7165b50f9d
Branch: yt
User: atmyers
Date: 2016-05-12 20:08:34+00:00
Summary: fixing a bug in unstructured mesh slicing when ax=1.
Affected #: 1 file
diff -r 6068831a37f4d92db1f3788c867e56efb286e587 -r 5d7165b50f9d014a67c32febac9e7a79e40a29f2 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -76,10 +76,15 @@
field_data = ad[field]
buff_size = size[0:dimension] + (1,) + size[dimension:]
+ ax = data_source.axis
+ xax = self.x_axis[ax]
+ yax = self.y_axis[ax]
c = np.float64(data_source.center[dimension].d)
- bounds.insert(2*dimension, c)
- bounds.insert(2*dimension, c)
- bounds = np.reshape(bounds, (3, 2))
+
+ extents = np.zeros((3, 2))
+ extents[ax] = np.array([c, c])
+ extents[xax] = bounds[0:2]
+ extents[yax] = bounds[2:4]
# if this is an element field, promote to 2D here
if len(field_data.shape) == 1:
@@ -101,13 +106,10 @@
img = pixelize_element_mesh(coords,
indices,
- buff_size, field_data, bounds,
+ buff_size, field_data, extents,
index_offset=offset)
# re-order the array and squeeze out the dummy dim
- ax = data_source.axis
- xax = self.x_axis[ax]
- yax = self.y_axis[ax]
return np.squeeze(np.transpose(img, (yax, xax, ax)))
elif dimension < 3:
https://bitbucket.org/yt_analysis/yt/commits/7abc3640e6c7/
Changeset: 7abc3640e6c7
Branch: yt
User: atmyers
Date: 2016-05-13 01:14:07+00:00
Summary: abort if triangle is parallel to plane
Affected #: 1 file
diff -r 5d7165b50f9d014a67c32febac9e7a79e40a29f2 -r 7abc3640e6c752d64239c9bad510208289dcbcf6 yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -18,8 +18,9 @@
cimport cython
from libc.stdlib cimport malloc, free
from yt.utilities.lib.fp_utils cimport fclip, i64clip
-from libc.math cimport copysign
+from libc.math cimport copysign, fabs
from yt.utilities.exceptions import YTDomainOverflow
+from yt.utilities.lib.vec3_ops cimport subtract, cross, dot, L2_norm
cdef extern from "math.h":
double exp(double x) nogil
@@ -501,6 +502,11 @@
cdef np.float64_t p0[3]
cdef np.float64_t p1[3]
cdef np.float64_t p3[3]
+ cdef np.float64_t E0[3]
+ cdef np.float64_t E1[3]
+ cdef np.float64_t tri_norm[3]
+ cdef np.float64_t plane_norm[3]
+ cdef np.float64_t dp
cdef int i, j, k, count, ntri, nlines
nlines = 0
ntri = triangles.shape[0]
@@ -510,6 +516,22 @@
first = last = points = NULL
for i in range(ntri):
count = 0
+
+ # skip if triangle is close to being parallel to plane
+ for j in range(3):
+ p0[j] = triangles[i, 0, j]
+ p1[j] = triangles[i, 1, j]
+ p3[j] = triangles[i, 2, j]
+ plane_norm[j] = 0.0
+ plane_norm[ax] = 1.0
+ subtract(p1, p0, E0)
+ subtract(p3, p0, E1)
+ cross(E0, E1, tri_norm)
+ dp = dot(tri_norm, plane_norm)
+ dp /= L2_norm(tri_norm)
+ if (fabs(dp) > 0.995):
+ continue
+
# Now for each line segment (01, 12, 20) we check to see how many cross
# the coordinate.
for j in range(3):
https://bitbucket.org/yt_analysis/yt/commits/3b1ff265bff3/
Changeset: 3b1ff265bff3
Branch: yt
User: atmyers
Date: 2016-05-13 06:47:25+00:00
Summary: use the triangle facet callback to perform mesh line annotation.
Affected #: 1 file
diff -r 7abc3640e6c752d64239c9bad510208289dcbcf6 -r 3b1ff265bff3ab218890c2286bc3e0074eb28e20 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -22,8 +22,6 @@
from distutils.version import LooseVersion
-from yt.config import \
- ytcfg
from yt.funcs import \
mylog, iterable
from yt.extern.six import add_metaclass
@@ -31,16 +29,37 @@
from yt.visualization.image_writer import apply_colormap
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.utilities.lib.pixelization_routines import \
- pixelize_element_mesh, pixelize_off_axis_cartesian, \
+ pixelize_off_axis_cartesian, \
pixelize_cartesian
from yt.analysis_modules.cosmological_observation.light_ray.light_ray import \
periodic_ray
from yt.utilities.lib.line_integral_convolution import \
line_integral_convolution_2d
-
+from yt.utilities.exceptions import YTElementTypeNotRecognized
from . import _MPL
+triangulate_hex = np.array([
+ [0, 2, 1], [0, 3, 2],
+ [4, 5, 6], [4, 6, 7],
+ [0, 1, 5], [0, 5, 4],
+ [1, 2, 6], [1, 6, 5],
+ [0, 7, 3], [0, 4, 7],
+ [3, 6, 2], [3, 7, 6]]
+)
+
+triangulate_tetra = np.array([
+ [0, 1, 3], [2, 3, 1],
+ [0, 3, 2], [0, 2, 1]]
+)
+
+triangulate_wedge = np.array([
+ [3, 0, 1], [4, 3, 1],
+ [2, 5, 4], [2, 4, 1],
+ [0, 3, 2], [2, 3, 5],
+ [3, 4, 5], [0, 2, 1]]
+)
+
callback_registry = {}
class RegisteredCallback(type):
@@ -1605,85 +1624,84 @@
class MeshLinesCallback(PlotCallback):
"""
- annotate_mesh_lines(thresh=0.01)
+ annotate_mesh_lines()
Adds the mesh lines to the plot.
- This is done by marking the pixels where the mapped coordinate
- is within some threshold distance of one of the element boundaries.
- If the mesh lines are too thick or too thin, try varying thresh.
-
Parameters
----------
- thresh : float
- The threshold distance, in mapped coordinates, within which the
- pixels will be marked as part of the element boundary.
- Default is 0.01.
+ plot_args: dict, optional
+ A dictionary of arguments that will be passed to matplotlib.
+
+ Example
+ -------
+
+ >>> import yt
+ >>> ds = yt.load("MOOSE_sample_data/out.e-s010")
+ >>> sl = yt.SlicePlot(ds, 'z', ('connect2', 'convected'))
+ >>> sl.annotate_mesh_lines(plot_args={'color':'black'})
"""
_type_name = "mesh_lines"
- def __init__(self, thresh=0.1, cmap=None):
+ def __init__(self, plot_args=None):
super(MeshLinesCallback, self).__init__()
- self.thresh = thresh
- if cmap is None:
- cmap = ytcfg.get("yt", "default_colormap")
- self.cmap = cmap
+ self.plot_args = plot_args
+
+ def promote_2d_to_3d(self, coords, indices, plot):
+ new_coords = np.zeros((2*coords.shape[0], 3))
+ new_connects = np.zeros((indices.shape[0], 2*indices.shape[1]),
+ dtype=np.int64)
+
+ new_coords[0:coords.shape[0],0:2] = coords
+ new_coords[0:coords.shape[0],2] = plot.ds.domain_left_edge[2]
+ new_coords[coords.shape[0]:,0:2] = coords
+ new_coords[coords.shape[0]:,2] = plot.ds.domain_right_edge[2]
+
+ new_connects[:,0:indices.shape[1]] = indices
+ new_connects[:,indices.shape[1]:] = indices + coords.shape[0]
+
+ return new_coords, new_connects
def __call__(self, plot):
+ index = plot.ds.index
ftype, fname = plot.field
mesh_id = int(ftype[-1]) - 1
- mesh = plot.ds.index.meshes[mesh_id]
+ mesh = index.meshes[mesh_id]
+ coords = mesh.connectivity_coords
+ indices = mesh.connectivity_indices - mesh._index_offset
- coords = mesh.connectivity_coords
- indices = mesh.connectivity_indices
+ num_verts = indices.shape[1]
+ num_dims = coords.shape[1]
- xx0, xx1 = plot._axes.get_xlim()
- yy0, yy1 = plot._axes.get_ylim()
+ if num_dims == 3 and (num_verts == 8 or
+ num_verts == 20 or
+ num_verts == 27):
+ tri_array = triangulate_hex
+ elif num_dims == 3 and num_verts == 4:
+ tri_array = triangulate_tetra
+ elif num_dims == 3 and num_verts == 6:
+ tri_array = triangulate_wedge
+ elif num_dims == 2 and num_verts == 3:
+ tri_array = triangulate_wedge
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
+ elif num_dims == 2 and num_verts == 4:
+ tri_array = triangulate_hex
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
+ else:
+ raise YTElementTypeNotRecognized(num_dims, num_verts)
- offset = mesh._index_offset
- ax = plot.data.axis
- xax = plot.data.ds.coordinates.x_axis[ax]
- yax = plot.data.ds.coordinates.y_axis[ax]
-
- size = plot.image._A.shape
- c = plot.data.center[ax]
- buff_size = size[0:ax] + (1,) + size[ax:]
-
- x0, x1 = plot.xlim
- y0, y1 = plot.ylim
- c = plot.data.center[ax]
- bounds = [x0, x1, y0, y1]
- bounds.insert(2*ax, c)
- bounds.insert(2*ax, c)
- bounds = np.reshape(bounds, (3, 2))
-
- # draw the mesh lines by marking where the mapped
- # coordinate is close to the domain edge. We do this by
- # calling the pixelizer with a dummy field and passing in
- # a non-negative thresh.
- dummy_field = np.zeros(indices.shape, dtype=np.float64)
- img = pixelize_element_mesh(coords, indices,
- buff_size,
- dummy_field,
- bounds,
- offset, self.thresh)
- img = np.squeeze(np.transpose(img, (yax, xax, ax)))
-
- # convert to RGBA
- size = plot.frb.buff_size
- image = np.zeros((size[0], size[1], 4), dtype=np.uint8)
- image[:, :, 0][img > 0.0] = 0
- image[:, :, 1][img > 0.0] = 0
- image[:, :, 2][img > 0.0] = 0
- image[:, :, 3][img > 0.0] = 255
-
- plot._axes.imshow(image, zorder=1,
- extent=[xx0, xx1, yy0, yy1],
- origin='lower', cmap=self.cmap,
- interpolation='nearest')
+ tri_indices = []
+ for elem in indices:
+ for tri in tri_array:
+ tri_indices.append(elem[tri])
+ tri_indices = np.array(tri_indices)
+ points = coords[tri_indices]
+
+ tfc = TriangleFacetsCallback(points, plot_args=self.plot_args)
+ return tfc(plot)
class TriangleFacetsCallback(PlotCallback):
https://bitbucket.org/yt_analysis/yt/commits/f9083b0b262c/
Changeset: f9083b0b262c
Branch: yt
User: atmyers
Date: 2016-05-13 06:48:32+00:00
Summary: remove some now un-needed code from pixelize_element_mesh.
Affected #: 1 file
diff -r 3b1ff265bff3ab218890c2286bc3e0074eb28e20 -r f9083b0b262c19c2ebb1231509f043ac238964a9 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -18,7 +18,9 @@
cimport cython
cimport libc.math as math
from yt.utilities.lib.fp_utils cimport fmin, fmax, i64min, i64max, imin, imax, fabs
-from yt.utilities.exceptions import YTPixelizeError, \
+from yt.utilities.exceptions import \
+ YTPixelizeError, \
+ YTException, \
YTElementTypeNotRecognized
from libc.stdlib cimport malloc, free
from vec3_ops cimport dot, cross, subtract
@@ -544,8 +546,7 @@
buff_size,
np.ndarray[np.float64_t, ndim=2] field,
extents,
- int index_offset = 0,
- double thresh = -1.0):
+ int index_offset = 0):
cdef np.ndarray[np.float64_t, ndim=3] img
img = np.zeros(buff_size, dtype="float64")
# Two steps:
@@ -592,27 +593,11 @@
# if we are in 2D land, the 1 cell thick dimension had better be 'z'
if ndim == 2:
- assert(buff_size[2] == 1)
+ try:
+ assert(buff_size[2] == 1)
+ except AssertionError:
+ raise YTException("2D datasets can only be sliced in the 'z' direction.")
- ax = -1
- for i in range(3):
- if buff_size[i] == 1:
- ax = i
- if ax == -1:
- raise RuntimeError
- xax = yax = -1
- if ax == 0:
- xax = 1
- yax = 2
- elif ax == 1:
- xax = 1
- yax = 2
- elif ax == 2:
- xax = 0
- yax = 1
- if xax == -1 or yax == -1:
- raise RuntimeError
-
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
@@ -673,20 +658,11 @@
sampler.map_real_to_unit(mapped_coord, vertices, ppoint)
if not sampler.check_inside(mapped_coord):
continue
- # if thresh is negative, we do the normal sample operation
- if thresh < 0.0:
- if (num_field_vals == 1):
- img[pi, pj, pk] = field_vals[0]
- else:
- img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
- field_vals)
+ if (num_field_vals == 1):
+ img[pi, pj, pk] = field_vals[0]
else:
- # otherwise, we draw the element boundaries
- if sampler.check_near_edge(mapped_coord, thresh, xax):
- img[pi, pj, pk] = 1.0
- elif sampler.check_near_edge(mapped_coord, thresh, yax):
- img[pi, pj, pk] = 1.0
-
+ img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
+ field_vals)
free(vertices)
free(field_vals)
return img
https://bitbucket.org/yt_analysis/yt/commits/242b0793fc56/
Changeset: 242b0793fc56
Branch: yt
User: atmyers
Date: 2016-05-13 06:49:17+00:00
Summary: remove the check_near_edge methods from the element samplers.
Affected #: 2 files
diff -r f9083b0b262c19c2ebb1231509f043ac238964a9 -r 242b0793fc566d2628f26ecce63c0519cbac31c7 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -30,11 +30,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -52,11 +47,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef class P1Sampler3D(ElementSampler):
@@ -72,11 +62,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -169,11 +154,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -191,11 +171,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -214,11 +189,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -250,8 +220,3 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
-
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
diff -r f9083b0b262c19c2ebb1231509f043ac238964a9 -r 242b0793fc566d2628f26ecce63c0519cbac31c7 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -86,15 +86,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- pass
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
pass
@@ -182,19 +173,6 @@
return 0
return 1
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
cdef class P1Sampler3D(ElementSampler):
'''
@@ -261,19 +239,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double u, v, w
cdef double thresh = 2.0e-2
@@ -413,18 +378,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
@@ -568,19 +521,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -781,22 +721,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (direction == 2 and (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance)):
- return 1
- elif (direction == 1 and (fabs(mapped_coord[direction]) < tolerance)):
- return 1
- elif (direction == 0 and (fabs(mapped_coord[0]) < tolerance or
- fabs(mapped_coord[0] + mapped_coord[1] - 1.0) < tolerance)):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double r, s, t
cdef double thresh = 5.0e-2
@@ -969,20 +893,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
https://bitbucket.org/yt_analysis/yt/commits/909926cb5d63/
Changeset: 909926cb5d63
Branch: yt
User: atmyers
Date: 2016-05-13 06:49:46+00:00
Summary: update the doc example of mesh line annotation on slices.
Affected #: 1 file
diff -r 242b0793fc566d2628f26ecce63c0519cbac31c7 -r 909926cb5d633b7e0969a1814d782a4a5ce912cd doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -406,16 +406,14 @@
import yt
ds = yt.load('MOOSE_sample_data/out.e-s010')
sl = yt.SlicePlot(ds, 'z', ('connect1', 'diffused'))
- sl.annotate_mesh_lines(thresh=0.1)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
sl.zoom(0.75)
sl.save()
-This annotation is performed by marking the pixels where the mapped coordinate is close
-to the element boundary. What counts as 'close' (in the mapped coordinate system) is
-determined by the ``thresh`` parameter, which can be varied to make the lines thicker or
-thinner.
+The ``plot_args`` parameter is a dictionary of keyword arguments that will be passed
+to matplotlib. It can be used to control the mesh line color, thickness, etc...
-The above example all involve 8-node hexahedral mesh elements. Here is another example from
+The above examples all involve 8-node hexahedral mesh elements. Here is another example from
a dataset that uses 6-node wedge elements:
.. python-script::
https://bitbucket.org/yt_analysis/yt/commits/eb9e446f1319/
Changeset: eb9e446f1319
Branch: yt
User: atmyers
Date: 2016-05-13 06:50:11+00:00
Summary: pixelization routines does need to be C++
Affected #: 1 file
diff -r 909926cb5d633b7e0969a1814d782a4a5ce912cd -r eb9e446f13191994fe3b92b6bfea64315504ebbe setup.py
--- a/setup.py
+++ b/setup.py
@@ -200,7 +200,6 @@
["yt/utilities/lib/pixelization_routines.pyx",
"yt/utilities/lib/pixelization_constants.c"],
include_dirs=["yt/utilities/lib/"],
- language="c++",
libraries=std_libs, depends=["yt/utilities/lib/fp_utils.pxd",
"yt/utilities/lib/pixelization_constants.h",
"yt/utilities/lib/element_mappings.pxd"]),
https://bitbucket.org/yt_analysis/yt/commits/3bd9cbb7594b/
Changeset: 3bd9cbb7594b
Branch: yt
User: atmyers
Date: 2016-05-23 16:50:44+00:00
Summary: merging with tip
Affected #: 41 files
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 .hgchurn
--- a/.hgchurn
+++ b/.hgchurn
@@ -1,6 +1,6 @@
stephenskory at yahoo.com = s at skory.us
"Stephen Skory stephenskory at yahoo.com" = s at skory.us
-yuan at astro.columbia.edu = bear0980 at gmail.com
+bear0980 at gmail.com = yuan at astro.columbia.edu
juxtaposicion at gmail.com = cemoody at ucsc.edu
chummels at gmail.com = chummels at astro.columbia.edu
jwise at astro.princeton.edu = jwise at physics.gatech.edu
@@ -19,7 +19,6 @@
sername=kayleanelson = kaylea.nelson at yale.edu
kayleanelson = kaylea.nelson at yale.edu
jcforbes at ucsc.edu = jforbes at ucolick.org
-ngoldbau at ucsc.edu = goldbaum at ucolick.org
biondo at wisc.edu = Biondo at wisc.edu
samgeen at googlemail.com = samgeen at gmail.com
fbogert = fbogert at ucsc.edu
@@ -39,4 +38,12 @@
jnaiman at ucolick.org = jnaiman
migueld.deval = miguel at archlinux.net
slevy at ncsa.illinois.edu = salevy at illinois.edu
-malzraa at gmail.com = kellerbw at mcmaster.ca
\ No newline at end of file
+malzraa at gmail.com = kellerbw at mcmaster.ca
+None = convert-repo
+dfenn = df11c at my.fsu.edu
+langmm = langmm.astro at gmail.com
+jmt354 = jmtomlinson95 at gmail.com
+desika = dnarayan at haverford.edu
+Ben Thompson = bthompson2090 at gmail.com
+goldbaum at ucolick.org = ngoldbau at illinois.edu
+ngoldbau at ucsc.edu = ngoldbau at illinois.edu
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 CREDITS
--- a/CREDITS
+++ b/CREDITS
@@ -10,6 +10,7 @@
Alex Bogert (fbogert at ucsc.edu)
André-Patrick Bubel (code at andre-bubel.de)
Pengfei Chen (madcpf at gmail.com)
+ Yi-Hao Chen (yihaochentw at gmail.com)
David Collins (dcollins4096 at gmail.com)
Brian Crosby (crosby.bd at gmail.com)
Andrew Cunningham (ajcunn at gmail.com)
@@ -25,10 +26,12 @@
William Gray (graywilliamj at gmail.com)
Markus Haider (markus.haider at uibk.ac.at)
Eric Hallman (hallman13 at gmail.com)
+ David Hannasch (David.A.Hannasch at gmail.com)
Cameron Hummels (chummels at gmail.com)
Anni Järvenpää (anni.jarvenpaa at gmail.com)
Allyson Julian (astrohckr at gmail.com)
Christian Karch (chiffre at posteo.de)
+ Maximilian Katz (maximilian.katz at stonybrook.edu)
Ben W. Keller (kellerbw at mcmaster.ca)
Ji-hoon Kim (me at jihoonkim.org)
Steffen Klemer (sklemer at phys.uni-goettingen.de)
@@ -60,6 +63,7 @@
Anna Rosen (rosen at ucolick.org)
Chuck Rozhon (rozhon2 at illinois.edu)
Douglas Rudd (drudd at uchicago.edu)
+ Hsi-Yu Schive (hyschive at gmail.com)
Anthony Scopatz (scopatz at gmail.com)
Noel Scudder (noel.scudder at stonybrook.edu)
Pat Shriwise (shriwise at wisc.edu)
@@ -75,6 +79,7 @@
Elizabeth Tasker (tasker at astro1.sci.hokudai.ac.jp)
Benjamin Thompson (bthompson2090 at gmail.com)
Robert Thompson (rthompsonj at gmail.com)
+ Joseph Tomlinson (jmtomlinson95 at gmail.com)
Stephanie Tonnesen (stonnes at gmail.com)
Matthew Turk (matthewturk at gmail.com)
Rich Wagner (rwagner at physics.ucsd.edu)
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -775,32 +775,38 @@
----------
FLASH HDF5 data is *mostly* supported and cared for by John ZuHone. To load a
-FLASH dataset, you can use the ``yt.load`` command and provide it the file name of a plot file or checkpoint file, but particle
-files are not currently directly loadable by themselves, due to the fact that
-they typically lack grid information. For instance, if you were in a directory
-with the following files:
+FLASH dataset, you can use the ``yt.load`` command and provide it the file name of
+a plot file, checkpoint file, or particle file. Particle files require special handling
+depending on the situation, the main issue being that they typically lack grid information.
+The first case is when you have a plotfile and a particle file that you would like to
+load together. In the simplest case, this occurs automatically. For instance, if you
+were in a directory with the following files:
.. code-block:: none
- cosmoSim_coolhdf5_chk_0026
+ radio_halo_1kpc_hdf5_plt_cnt_0100 # plotfile
+ radio_halo_1kpc_hdf5_part_0100 # particle file
-You would feed it the filename ``cosmoSim_coolhdf5_chk_0026``:
+where the plotfile and the particle file were created at the same time (therefore having
+particle data consistent with the grid structure of the former). Notice also that the
+prefix ``"radio_halo_1kpc_"`` and the file number ``100`` are the same. In this special case,
+the particle file will be loaded automatically when ``yt.load`` is called on the plotfile.
+This also works when loading a number of files in a time series.
-.. code-block:: python
-
- import yt
- ds = yt.load("cosmoSim_coolhdf5_chk_0026")
-
-If you have a FLASH particle file that was created at the same time as
-a plotfile or checkpoint file (therefore having particle data
-consistent with the grid structure of the latter), its data may be loaded with the
-``particle_filename`` optional argument:
+If the two files do not have the same prefix and number, but they nevertheless have the same
+grid structure and are at the same simulation time, the particle data may be loaded with the
+``particle_filename`` optional argument to ``yt.load``:
.. code-block:: python
import yt
ds = yt.load("radio_halo_1kpc_hdf5_plt_cnt_0100", particle_filename="radio_halo_1kpc_hdf5_part_0100")
+However, if you don't have a corresponding plotfile for a particle file, but would still
+like to load the particle data, you can still call ``yt.load`` on the file. However, the
+grid information will not be available, and the particle data will be loaded in a fashion
+similar to SPH data.
+
.. rubric:: Caveats
* Please be careful that the units are correctly utilized; yt assumes cgs by default, but conversion to
@@ -1313,6 +1319,31 @@
``bbox``
The bounding box for the particle positions.
+.. _loading-gizmo-data:
+
+Gizmo Data
+----------
+
+Gizmo datasets, including FIRE outputs, can be loaded into yt in the usual
+manner. Like other SPH data formats, yt loads Gizmo data as particle fields
+and then uses smoothing kernels to deposit those fields to an underlying
+grid structure as spatial fields as described in :ref:`loading-gadget-data`.
+To load Gizmo datasets using the standard HDF5 output format::
+
+.. code-block:: python
+
+ import yt
+ ds = yt.load("snapshot_600.hdf5")
+
+Because the Gizmo output format is similar to the Gadget format, yt
+may load Gizmo datasets as Gadget depending on the circumstances, but this
+should not pose a problem in most situations. FIRE outputs will be loaded
+accordingly due to the number of metallicity fields found (11 or 17).
+
+For Gizmo outputs written as raw binary outputs, you may have to specify
+a bounding box, field specification, and units as are done for standard
+Gadget outputs. See :ref:`loading-gadget-data` for more information.
+
.. _loading-pyne-data:
Halo Catalog Data
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 doc/source/reference/code_support.rst
--- a/doc/source/reference/code_support.rst
+++ b/doc/source/reference/code_support.rst
@@ -28,16 +28,18 @@
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| Enzo | Y | Y | Y | Y | Y | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
+| FITS | Y | N/A | Y | Y | Y | Y | Y | Full |
++-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| FLASH | Y | Y | Y | Y | Y | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
-| FITS | Y | N/A | Y | Y | Y | Y | Y | Full |
-+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| Gadget | Y | Y | Y | Y | Y [#f2]_ | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| GAMER | Y | N | Y | Y | Y | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| Gasoline | Y | Y | Y | Y | Y [#f2]_ | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
+| Gizmo | Y | Y | Y | Y | Y [#f2]_ | Y | Y | Full |
++-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| Grid Data Format (GDF)| Y | N/A | Y | Y | Y | Y | Y | Full |
+-----------------------+------------+-----------+------------+-------+----------+----------+------------+----------+
| Maestro | Y [#f1]_ | N | Y | Y | Y | Y | N | Partial |
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 setup.py
--- a/setup.py
+++ b/setup.py
@@ -451,7 +451,7 @@
license="BSD",
zip_safe=False,
scripts=["scripts/iyt"],
- data_files=MAPSERVER_FILES + SHADERS_FILES,
+ data_files=MAPSERVER_FILES + [(SHADERS_DIR, SHADERS_FILES)],
ext_modules=cython_extensions + extensions
)
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -14,7 +14,7 @@
local_fits_000:
- yt/frontends/fits/tests/test_outputs.py
- local_flash_000:
+ local_flash_001:
- yt/frontends/flash/tests/test_outputs.py
local_gadget_000:
@@ -26,6 +26,9 @@
local_gdf_000:
- yt/frontends/gdf/tests/test_outputs.py
+ local_gizmo_000:
+ - yt/frontends/gizmo/tests/test_outputs.py
+
local_halos_000:
- yt/analysis_modules/halo_analysis/tests/test_halo_finders.py # [py2]
- yt/analysis_modules/halo_finding/tests/test_rockstar.py # [py2]
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
--- a/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
+++ b/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
@@ -35,10 +35,12 @@
class LightRay(CosmologySplice):
"""
- A LightRay object is a one-dimensional object representing the trajectory
- of a ray of light as it passes through one or more datasets (simple and
- compound rays respectively). One can sample any of the fields intersected
- by the LightRay object as it passed through the dataset(s).
+ A 1D object representing the path of a light ray passing through a
+ simulation. LightRays can be either simple, where they pass through a
+ single dataset, or compound, where they pass through consecutive
+ datasets from the same cosmological simulation. One can sample any of
+ the fields intersected by the LightRay object as it passed through
+ the dataset(s).
For compound rays, the LightRay stacks together multiple datasets in a time
series in order to approximate a LightRay's path through a volume
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/data_objects/octree_subset.py
--- a/yt/data_objects/octree_subset.py
+++ b/yt/data_objects/octree_subset.py
@@ -433,57 +433,62 @@
self.base_region = base_region
self.base_selector = base_region.selector
+class OctreeSubsetBlockSlicePosition(object):
+ def __init__(self, ind, block_slice):
+ self.ind = ind
+ self.block_slice = block_slice
+ nz = self.block_slice.octree_subset.nz
+ self.ActiveDimensions = np.array([nz,nz,nz], dtype="int64")
+
+ def __getitem__(self, key):
+ bs = self.block_slice
+ rv = bs.octree_subset[key][:,:,:,self.ind]
+ if bs.octree_subset._block_reorder:
+ rv = rv.copy(order=bs.octree_subset._block_reorder)
+ return rv
+
+ @property
+ def id(self):
+ return self.ind
+
+ @property
+ def Level(self):
+ return self.block_slice._ires[0,0,0,self.ind]
+
+ @property
+ def LeftEdge(self):
+ LE = (self.block_slice._fcoords[0,0,0,self.ind,:]
+ - self.block_slice._fwidth[0,0,0,self.ind,:]*0.5)
+ return LE
+
+ @property
+ def RightEdge(self):
+ RE = (self.block_slice._fcoords[-1,-1,-1,self.ind,:]
+ + self.block_slice._fwidth[-1,-1,-1,self.ind,:]*0.5)
+ return RE
+
+ @property
+ def dds(self):
+ return self.block_slice._fwidth[0,0,0,self.ind,:]
+
+ def clear_data(self):
+ pass
+
+ def get_vertex_centered_data(self, *args, **kwargs):
+ raise NotImplementedError
+
+
class OctreeSubsetBlockSlice(object):
def __init__(self, octree_subset):
- self.ind = None
self.octree_subset = octree_subset
# Cache some attributes
- nz = octree_subset.nz
- self.ActiveDimensions = np.array([nz,nz,nz], dtype="int64")
for attr in ["ires", "icoords", "fcoords", "fwidth"]:
v = getattr(octree_subset, attr)
setattr(self, "_%s" % attr, octree_subset._reshape_vals(v))
def __iter__(self):
for i in range(self._ires.shape[-1]):
- self.ind = i
- yield i, self
-
- def clear_data(self):
- pass
-
- def __getitem__(self, key):
- rv = self.octree_subset[key][:,:,:,self.ind]
- if self.octree_subset._block_reorder:
- rv = rv.copy(order=self.octree_subset._block_reorder)
- return rv
-
- def get_vertex_centered_data(self, *args, **kwargs):
- raise NotImplementedError
-
- @property
- def id(self):
- return np.random.randint(1)
-
- @property
- def Level(self):
- return self._ires[0,0,0,self.ind]
-
- @property
- def LeftEdge(self):
- LE = (self._fcoords[0,0,0,self.ind,:]
- - self._fwidth[0,0,0,self.ind,:]*0.5)
- return LE
-
- @property
- def RightEdge(self):
- RE = (self._fcoords[-1,-1,-1,self.ind,:]
- + self._fwidth[-1,-1,-1,self.ind,:]*0.5)
- return RE
-
- @property
- def dds(self):
- return self._fwidth[0,0,0,self.ind,:]
+ yield i, OctreeSubsetBlockSlicePosition(i, self)
class YTPositionArray(YTArray):
@property
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/data_objects/unstructured_mesh.py
--- a/yt/data_objects/unstructured_mesh.py
+++ b/yt/data_objects/unstructured_mesh.py
@@ -225,3 +225,12 @@
else:
self._last_count = mask.sum()
return mask
+
+ def select(self, selector, source, dest, offset):
+ mask = self._get_selector_mask(selector)
+ count = self.count(selector)
+ if count == 0: return 0
+ # Note: this likely will not work with vector fields.
+ dest[offset:offset+count] = source.flat[mask]
+ return count
+
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/api.py
--- a/yt/frontends/api.py
+++ b/yt/frontends/api.py
@@ -31,6 +31,7 @@
'gadget_fof',
'gamer',
'gdf',
+ 'gizmo',
'halo_catalog',
'http_stream',
'moab',
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -444,6 +444,16 @@
if "AppendActiveParticleType" in self.dataset.parameters:
ap_fields = self._detect_active_particle_fields()
field_list = list(set(field_list).union(ap_fields))
+ if not any(f[0] == 'io' for f in field_list):
+ if 'io' in self.dataset.particle_types_raw:
+ ptypes_raw = list(self.dataset.particle_types_raw)
+ ptypes_raw.remove('io')
+ self.dataset.particle_types_raw = tuple(ptypes_raw)
+
+ if 'io' in self.dataset.particle_types:
+ ptypes = list(self.dataset.particle_types)
+ ptypes.remove('io')
+ self.dataset.particle_types = tuple(ptypes)
ptypes = self.dataset.particle_types
ptypes_raw = self.dataset.particle_types_raw
else:
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/enzo/tests/test_outputs.py
--- a/yt/frontends/enzo/tests/test_outputs.py
+++ b/yt/frontends/enzo/tests/test_outputs.py
@@ -126,7 +126,9 @@
def test_active_particle_datasets():
two_sph = data_dir_load(two_sphere_test)
assert 'AccretingParticle' in two_sph.particle_types_raw
- assert_equal(len(two_sph.particle_unions), 0)
+ assert 'io' not in two_sph.particle_types_raw
+ assert 'all' in two_sph.particle_types
+ assert_equal(len(two_sph.particle_unions), 1)
pfields = ['GridID', 'creation_time', 'dynamical_time',
'identifier', 'level', 'metallicity', 'particle_mass']
pfields += ['particle_position_%s' % d for d in 'xyz']
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/exodus_ii/data_structures.py
--- a/yt/frontends/exodus_ii/data_structures.py
+++ b/yt/frontends/exodus_ii/data_structures.py
@@ -154,7 +154,8 @@
units_override=units_override)
self.index_filename = filename
self.storage_filename = storage_filename
- self.default_field = ("connect1", "diffused")
+ self.default_field = [f for f in self.field_list
+ if f[0] == 'connect1'][-1]
def _set_code_unit_attributes(self):
# This is where quantities are created that represent the various
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/flash/api.py
--- a/yt/frontends/flash/api.py
+++ b/yt/frontends/flash/api.py
@@ -16,12 +16,14 @@
from .data_structures import \
FLASHGrid, \
FLASHHierarchy, \
- FLASHDataset
+ FLASHDataset, \
+ FLASHParticleDataset
from .fields import \
FLASHFieldInfo
from .io import \
- IOHandlerFLASH
+ IOHandlerFLASH, \
+ IOHandlerFLASHParticle
from . import tests
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/flash/data_structures.py
--- a/yt/frontends/flash/data_structures.py
+++ b/yt/frontends/flash/data_structures.py
@@ -18,13 +18,15 @@
import numpy as np
import weakref
-from yt.funcs import mylog
from yt.data_objects.grid_patch import \
AMRGridPatch
+from yt.data_objects.static_output import \
+ Dataset, ParticleFile
+from yt.funcs import mylog
from yt.geometry.grid_geometry_handler import \
GridIndex
-from yt.data_objects.static_output import \
- Dataset
+from yt.geometry.particle_geometry_handler import \
+ ParticleIndex
from yt.utilities.file_handler import \
HDF5FileHandler
from yt.utilities.physical_ratios import cm_per_mpc
@@ -251,8 +253,6 @@
self.time_unit = self.quan(1.0, "s")
self.velocity_unit = self.quan(1.0, "cm/s")
self.temperature_unit = self.quan(temperature_factor, "K")
- # Still need to deal with:
- #self.conversion_factors['temp'] = (1.0 + self.current_redshift)**-2.0
self.unit_registry.modify("code_magnetic", self.magnetic_unit)
def set_code_units(self):
@@ -437,3 +437,52 @@
def close(self):
self._handle.close()
+
+class FLASHParticleFile(ParticleFile):
+ pass
+
+class FLASHParticleDataset(FLASHDataset):
+ _index_class = ParticleIndex
+ over_refine_factor = 1
+ filter_bbox = False
+ _file_class = FLASHParticleFile
+
+ def __init__(self, filename, dataset_type='flash_particle_hdf5',
+ storage_filename = None,
+ units_override = None,
+ n_ref = 64, unit_system = "cgs"):
+
+ if self._handle is not None: return
+ self._handle = HDF5FileHandler(filename)
+ self.n_ref = n_ref
+ self.refine_by = 2
+ Dataset.__init__(self, filename, dataset_type, units_override=units_override,
+ unit_system=unit_system)
+ self.storage_filename = storage_filename
+
+ def _parse_parameter_file(self):
+ # Let the superclass do all the work but then
+ # fix the domain dimensions
+ super(FLASHParticleDataset, self)._parse_parameter_file()
+ nz = 1 << self.over_refine_factor
+ self.domain_dimensions = np.zeros(3, "int32")
+ self.domain_dimensions[:self.dimensionality] = nz
+ self.filename_template = self.parameter_filename
+ self.file_count = 1
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = HDF5FileHandler(args[0])
+ if "bounding box" not in fileh["/"].keys() \
+ and "localnp" in fileh["/"].keys():
+ return True
+ except IOError:
+ pass
+ return False
+
+ @classmethod
+ def _guess_candidates(cls, base, directories, files):
+ candidates = [_ for _ in files if "_hdf5_part_" in _]
+ # Typically, Flash won't have nested outputs.
+ return candidates, (len(candidates) == 0)
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/flash/io.py
--- a/yt/frontends/flash/io.py
+++ b/yt/frontends/flash/io.py
@@ -20,6 +20,8 @@
BaseIOHandler
from yt.utilities.logger import ytLogger as mylog
from yt.geometry.selection_routines import AlwaysSelector
+from yt.utilities.lib.geometry_utils import \
+ compute_morton
# http://stackoverflow.com/questions/2361945/detecting-consecutive-integers-in-a-list
def particle_sequences(grids):
@@ -34,6 +36,16 @@
seq = list(v[1] for v in g)
yield seq
+def determine_particle_fields(handle):
+ try:
+ particle_fields = [s[0].decode("ascii","ignore").strip()
+ for s in handle["/particle names"][:]]
+ _particle_fields = dict([("particle_" + s, i) for i, s in
+ enumerate(particle_fields)])
+ except KeyError:
+ _particle_fields = {}
+ return _particle_fields
+
class IOHandlerFLASH(BaseIOHandler):
_particle_reader = False
_dataset_type = "flash_hdf5"
@@ -43,15 +55,7 @@
# Now we cache the particle fields
self._handle = ds._handle
self._particle_handle = ds._particle_handle
-
- try :
- particle_fields = [s[0].decode("ascii","ignore").strip()
- for s in
- self._particle_handle["/particle names"][:]]
- self._particle_fields = dict([("particle_" + s, i) for i, s in
- enumerate(particle_fields)])
- except KeyError:
- self._particle_fields = {}
+ self._particle_fields = determine_particle_fields(self._particle_handle)
def _read_particles(self, fields_to_read, type, args, grid_list,
count_list, conv_factors):
@@ -154,3 +158,96 @@
rv[g.id][field] = np.asarray(data[...,i], "=f8")
return rv
+class IOHandlerFLASHParticle(BaseIOHandler):
+ _particle_reader = True
+ _dataset_type = "flash_particle_hdf5"
+
+ def __init__(self, ds):
+ super(IOHandlerFLASHParticle, self).__init__(ds)
+ # Now we cache the particle fields
+ self._handle = ds._handle
+ self._particle_fields = determine_particle_fields(self._handle)
+ self._position_fields = [self._particle_fields["particle_pos%s" % ax]
+ for ax in 'xyz']
+ self._chunksize = 32**3
+
+ def _read_fluid_selection(self, chunks, selector, fields, size):
+ raise NotImplementedError
+
+ def _read_particle_coords(self, chunks, ptf):
+ chunks = list(chunks)
+ data_files = set([])
+ assert(len(ptf) == 1)
+ for chunk in chunks:
+ for obj in chunk.objs:
+ data_files.update(obj.data_files)
+ px, py, pz = self._position_fields
+ p_fields = self._handle["/tracer particles"]
+ assert(len(data_files) == 1)
+ for data_file in sorted(data_files):
+ pcount = self._count_particles(data_file)["io"]
+ for ptype, field_list in sorted(ptf.items()):
+ total = 0
+ while total < pcount:
+ count = min(self._chunksize, pcount - total)
+ x = np.asarray(p_fields[total:total+count, px], dtype="=f8")
+ y = np.asarray(p_fields[total:total+count, py], dtype="=f8")
+ z = np.asarray(p_fields[total:total+count, pz], dtype="=f8")
+ total += count
+ yield ptype, (x, y, z)
+
+ def _read_particle_fields(self, chunks, ptf, selector):
+ chunks = list(chunks)
+ data_files = set([])
+ assert(len(ptf) == 1)
+ for chunk in chunks:
+ for obj in chunk.objs:
+ data_files.update(obj.data_files)
+ px, py, pz = self._position_fields
+ p_fields = self._handle["/tracer particles"]
+ assert(len(data_files) == 1)
+ for data_file in sorted(data_files):
+ pcount = self._count_particles(data_file)["io"]
+ for ptype, field_list in sorted(ptf.items()):
+ total = 0
+ while total < pcount:
+ count = min(self._chunksize, pcount - total)
+ x = np.asarray(p_fields[total:total+count, px], dtype="=f8")
+ y = np.asarray(p_fields[total:total+count, py], dtype="=f8")
+ z = np.asarray(p_fields[total:total+count, pz], dtype="=f8")
+ total += count
+ mask = selector.select_points(x, y, z, 0.0)
+ del x, y, z
+ if mask is None: continue
+ for field in field_list:
+ fi = self._particle_fields[field]
+ data = p_fields[total-count:total, fi]
+ yield (ptype, field), data[mask]
+
+ def _initialize_index(self, data_file, regions):
+ p_fields = self._handle["/tracer particles"]
+ px, py, pz = self._position_fields
+ pcount = self._count_particles(data_file)["io"]
+ morton = np.empty(pcount, dtype='uint64')
+ ind = 0
+ while ind < pcount:
+ npart = min(self._chunksize, pcount - ind)
+ pos = np.empty((npart, 3), dtype="=f8")
+ pos[:,0] = p_fields[ind:ind+npart, px]
+ pos[:,1] = p_fields[ind:ind+npart, py]
+ pos[:,2] = p_fields[ind:ind+npart, pz]
+ regions.add_data_file(pos, data_file.file_id)
+ morton[ind:ind+npart] = \
+ compute_morton(pos[:,0], pos[:,1], pos[:,2],
+ data_file.ds.domain_left_edge,
+ data_file.ds.domain_right_edge)
+ ind += self._chunksize
+ return morton
+
+ def _count_particles(self, data_file):
+ pcount = {"io": self._handle["/localnp"][:].sum()}
+ return pcount
+
+ def _identify_fields(self, data_file):
+ fields = [("io", field) for field in self._particle_fields]
+ return fields, {}
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/flash/tests/test_outputs.py
--- a/yt/frontends/flash/tests/test_outputs.py
+++ b/yt/frontends/flash/tests/test_outputs.py
@@ -20,8 +20,11 @@
from yt.utilities.answer_testing.framework import \
requires_ds, \
small_patch_amr, \
- data_dir_load
-from yt.frontends.flash.api import FLASHDataset
+ data_dir_load, \
+ sph_answer
+from yt.frontends.flash.api import FLASHDataset, \
+ FLASHParticleDataset
+from collections import OrderedDict
_fields = ("temperature", "density", "velocity_magnitude")
@@ -45,7 +48,6 @@
test_wind_tunnel.__name__ = test.description
yield test
-
@requires_file(wt)
def test_FLASHDataset():
assert isinstance(data_dir_load(wt), FLASHDataset)
@@ -54,3 +56,28 @@
def test_units_override():
for test in units_override_check(sloshing):
yield test
+
+fid_1to3_b1 = "fiducial_1to3_b1/fiducial_1to3_b1_hdf5_part_0080"
+
+fid_1to3_b1_fields = OrderedDict(
+ [
+ (("deposit", "all_density"), None),
+ (("deposit", "all_count"), None),
+ (("deposit", "all_cic"), None),
+ (("deposit", "all_cic_velocity_x"), ("deposit", "all_cic")),
+ (("deposit", "all_cic_velocity_y"), ("deposit", "all_cic")),
+ (("deposit", "all_cic_velocity_z"), ("deposit", "all_cic")),
+ ]
+)
+
+
+ at requires_file(fid_1to3_b1)
+def test_FLASHParticleDataset():
+ assert isinstance(data_dir_load(fid_1to3_b1), FLASHParticleDataset)
+
+ at requires_ds(fid_1to3_b1, big_data=True)
+def test_fid_1to3_b1():
+ ds = data_dir_load(fid_1to3_b1)
+ for test in sph_answer(ds, 'fiducial_1to3_b1_hdf5_part_0080', 6684119, fid_1to3_b1_fields):
+ test_fid_1to3_b1.__name__ = test.description
+ yield test
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/gizmo/api.py
--- /dev/null
+++ b/yt/frontends/gizmo/api.py
@@ -0,0 +1,21 @@
+"""
+API for Gizmo frontend.
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from .data_structures import \
+ GizmoDataset
+
+from .fields import \
+ GizmoFieldInfo
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/gizmo/data_structures.py
--- /dev/null
+++ b/yt/frontends/gizmo/data_structures.py
@@ -0,0 +1,44 @@
+"""
+Data structures for Gizmo frontend.
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from yt.utilities.on_demand_imports import _h5py as h5py
+
+from yt.frontends.gadget.data_structures import \
+ GadgetHDF5Dataset
+
+from .fields import \
+ GizmoFieldInfo
+
+class GizmoDataset(GadgetHDF5Dataset):
+ _field_info_class = GizmoFieldInfo
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ need_groups = ['Header']
+ veto_groups = ['FOF', 'Group', 'Subhalo']
+ valid = True
+ try:
+ fh = h5py.File(args[0], mode='r')
+ valid = all(ng in fh["/"] for ng in need_groups) and \
+ not any(vg in fh["/"] for vg in veto_groups)
+ dmetal = "/PartType0/Metallicity"
+ if dmetal not in fh or fh[dmetal].shape[1] not in (11, 17):
+ valid = False
+ fh.close()
+ except:
+ valid = False
+ pass
+ return valid
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/gizmo/fields.py
--- /dev/null
+++ b/yt/frontends/gizmo/fields.py
@@ -0,0 +1,116 @@
+"""
+Gizmo-specific fields
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from yt.fields.particle_fields import \
+ add_volume_weighted_smoothed_field
+from yt.fields.species_fields import \
+ add_species_field_by_density
+from yt.frontends.gadget.fields import \
+ GadgetFieldInfo
+from yt.frontends.sph.fields import \
+ SPHFieldInfo
+
+class GizmoFieldInfo(GadgetFieldInfo):
+ known_particle_fields = (
+ ("Mass", ("code_mass", ["particle_mass"], None)),
+ ("Masses", ("code_mass", ["particle_mass"], None)),
+ ("Coordinates", ("code_length", ["particle_position"], None)),
+ ("Velocity", ("code_velocity", ["particle_velocity"], None)),
+ ("Velocities", ("code_velocity", ["particle_velocity"], None)),
+ ("ParticleIDs", ("", ["particle_index"], None)),
+ ("InternalEnergy", ("code_velocity ** 2", ["thermal_energy"], None)),
+ ("SmoothingLength", ("code_length", ["smoothing_length"], None)),
+ ("Density", ("code_mass / code_length**3", ["density"], None)),
+ ("MaximumTemperature", ("K", [], None)),
+ ("Temperature", ("K", ["temperature"], None)),
+ ("Epsilon", ("code_length", [], None)),
+ ("Metals", ("code_metallicity", ["metallicity"], None)),
+ ("Metallicity", ("code_metallicity", ["metallicity"], None)),
+ ("Phi", ("code_length", [], None)),
+ ("StarFormationRate", ("Msun / yr", [], None)),
+ ("FormationTime", ("code_time", ["creation_time"], None)),
+ ("Metallicity_00", ("", ["metallicity"], None)),
+ ("Metallicity_01", ("", ["He_metallicity"], None)),
+ ("Metallicity_02", ("", ["C_metallicity"], None)),
+ ("Metallicity_03", ("", ["N_metallicity"], None)),
+ ("Metallicity_04", ("", ["O_metallicity"], None)),
+ ("Metallicity_05", ("", ["Ne_metallicity"], None)),
+ ("Metallicity_06", ("", ["Mg_metallicity"], None)),
+ ("Metallicity_07", ("", ["Si_metallicity"], None)),
+ ("Metallicity_08", ("", ["S_metallicity"], None)),
+ ("Metallicity_09", ("", ["Ca_metallicity"], None)),
+ ("Metallicity_10", ("", ["Fe_metallicity"], None)),
+ )
+
+ def __init__(self, *args, **kwargs):
+ super(SPHFieldInfo, self).__init__(*args, **kwargs)
+ if ("PartType0", "Metallicity_00") in self.field_list:
+ self.nuclei_names = ["He", "C", "N", "O", "Ne", "Mg", "Si", "S",
+ "Ca", "Fe"]
+
+ def setup_gas_particle_fields(self, ptype):
+ super(GizmoFieldInfo, self).setup_gas_particle_fields(ptype)
+ self.alias((ptype, "temperature"), (ptype, "Temperature"))
+
+ def _h_density(field, data):
+ x_H = 1.0 - data[(ptype, "He_metallicity")] - \
+ data[(ptype, "metallicity")]
+ return x_H * data[(ptype, "density")] * \
+ data[(ptype, "NeutralHydrogenAbundance")]
+
+ self.add_field(
+ (ptype, "H_density"),
+ function=_h_density,
+ particle_type=True,
+ units=self.ds.unit_system["density"])
+ add_species_field_by_density(self, ptype, "H", particle_type=True)
+ for suffix in ["density", "fraction", "mass", "number_density"]:
+ self.alias((ptype, "H_p0_%s" % suffix), (ptype, "H_%s" % suffix))
+
+ def _h_p1_density(field, data):
+ x_H = 1.0 - data[(ptype, "He_metallicity")] - \
+ data[(ptype, "metallicity")]
+ return x_H * data[(ptype, "density")] * \
+ (1.0 - data[(ptype, "NeutralHydrogenAbundance")])
+
+ self.add_field(
+ (ptype, "H_p1_density"),
+ function=_h_p1_density,
+ particle_type=True,
+ units=self.ds.unit_system["density"])
+ add_species_field_by_density(self, ptype, "H_p1", particle_type=True)
+
+ def _nuclei_mass_density_field(field, data):
+ species = field.name[1][:field.name[1].find("_")]
+ return data[ptype, "density"] * \
+ data[ptype, "%s_metallicity" % species]
+
+ num_neighbors = 64
+ for species in self.nuclei_names:
+ self.add_field(
+ (ptype, "%s_nuclei_mass_density" % species),
+ function=_nuclei_mass_density_field,
+ particle_type=True,
+ units=self.ds.unit_system["density"])
+
+ for suf in ["_nuclei_mass_density", "_metallicity"]:
+ field = "%s%s" % (species, suf)
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", field,
+ self, num_neighbors)
+
+ self.alias(("gas", field), fn[0])
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/gizmo/tests/test_outputs.py
--- /dev/null
+++ b/yt/frontends/gizmo/tests/test_outputs.py
@@ -0,0 +1,46 @@
+"""
+Gizmo frontend tests
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2015, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from collections import OrderedDict
+
+from yt.utilities.answer_testing.framework import \
+ data_dir_load, \
+ requires_ds, \
+ sph_answer
+from yt.frontends.gizmo.api import GizmoDataset
+
+FIRE_m12i = 'FIRE_M12i_ref11/snapshot_600.hdf5'
+
+# This maps from field names to weight field names to use for projections
+fields = OrderedDict(
+ [
+ (("gas", "density"), None),
+ (("gas", "temperature"), ('gas', 'density')),
+ (("gas", "metallicity"), ('gas', 'density')),
+ (("gas", "O_metallicity"), ('gas', 'density')),
+ (('gas', 'velocity_magnitude'), None),
+ (("deposit", "all_count"), None),
+ (("deposit", "all_cic"), None),
+ ]
+)
+
+ at requires_ds(FIRE_m12i)
+def test_GizmoDataset():
+ ds = data_dir_load(FIRE_m12i)
+ assert isinstance(ds, GizmoDataset)
+ for test in sph_answer(ds, 'snapshot_600', 4786950, fields):
+ test_GizmoDataset.__name__ = test.description
+ yield test
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -41,27 +41,9 @@
("Phi", ("code_length", [], None)),
("StarFormationRate", ("Msun / yr", [], None)),
("FormationTime", ("code_time", ["creation_time"], None)),
- # These are metallicity fields that get discovered for FIRE simulations
("Metallicity_00", ("", ["metallicity"], None)),
- ("Metallicity_01", ("", ["He_fraction"], None)),
- ("Metallicity_02", ("", ["C_fraction"], None)),
- ("Metallicity_03", ("", ["N_fraction"], None)),
- ("Metallicity_04", ("", ["O_fraction"], None)),
- ("Metallicity_05", ("", ["Ne_fraction"], None)),
- ("Metallicity_06", ("", ["Mg_fraction"], None)),
- ("Metallicity_07", ("", ["Si_fraction"], None)),
- ("Metallicity_08", ("", ["S_fraction"], None)),
- ("Metallicity_09", ("", ["Ca_fraction"], None)),
- ("Metallicity_10", ("", ["Fe_fraction"], None)),
)
- def __init__(self, *args, **kwargs):
- super(SPHFieldInfo, self).__init__(*args, **kwargs)
- # Special case for FIRE
- if ("PartType0", "Metallicity_00") in self.field_list:
- self.species_names += ["He", "C", "N", "O", "Ne", "Mg", "Si", "S",
- "Ca", "Fe"]
-
def setup_particle_fields(self, ptype, *args, **kwargs):
super(SPHFieldInfo, self).setup_particle_fields(ptype, *args, **kwargs)
setup_species_fields(self, ptype)
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -1826,5 +1826,7 @@
sds._node_fields = node_data[0].keys()
sds._elem_fields = elem_data[0].keys()
+ sds.default_field = [f for f in sds.field_list
+ if f[0] == 'connect1'][-1]
return sds
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/frontends/stream/tests/test_stream_hexahedral.py
--- a/yt/frontends/stream/tests/test_stream_hexahedral.py
+++ b/yt/frontends/stream/tests/test_stream_hexahedral.py
@@ -6,6 +6,7 @@
from yt.testing import \
assert_almost_equal, \
assert_equal
+from yt import SlicePlot
# Field information
def test_stream_hexahedral() :
@@ -30,7 +31,7 @@
cell_z[0] = 0.0
coords, conn = hexahedral_connectivity(cell_x, cell_y, cell_z)
- data = {'random_field': np.random.random(Nx*Ny*Nz)}
+ data = {'random_field': np.random.random((Nx, Ny, Nz))}
bbox = np.array([ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0] ])
ds = load_hexahedral_mesh(data, conn, coords, bbox=bbox)
dd = ds.all_data()
@@ -42,7 +43,7 @@
cell_y = np.linspace(0.0, 1.0, Ny+1)
cell_z = np.linspace(0.0, 1.0, Nz+1)
coords, conn = hexahedral_connectivity(cell_x, cell_y, cell_z)
- data = {'random_field': np.random.random(Nx*Ny*Nz)}
+ data = {'random_field': np.random.random((Nx, Ny, Nz))}
bbox = np.array([ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0] ])
ds = load_hexahedral_mesh(data, conn, coords, bbox=bbox)
dd = ds.all_data()
@@ -51,3 +52,7 @@
yield assert_almost_equal, dd["dx"].to_ndarray(), 1.0/Nx
yield assert_almost_equal, dd["dy"].to_ndarray(), 1.0/Ny
yield assert_almost_equal, dd["dz"].to_ndarray(), 1.0/Nz
+
+ s = SlicePlot(ds, "x", "random_field")
+ s._setup_plots()
+ s.frb["random_field"]
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/funcs.py
--- a/yt/funcs.py
+++ b/yt/funcs.py
@@ -318,7 +318,10 @@
ipshell(header = __header % dd,
local_ns = loc, global_ns = glo)
else:
- from IPython.config.loader import Config
+ try:
+ from traitlets.config.loader import Config
+ except ImportError:
+ from IPython.config.loader import Config
cfg = Config()
cfg.InteractiveShellEmbed.local_ns = loc
cfg.InteractiveShellEmbed.global_ns = glo
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/units/tests/test_ytarray.py
--- a/yt/units/tests/test_ytarray.py
+++ b/yt/units/tests/test_ytarray.py
@@ -112,6 +112,20 @@
yield assert_raises, YTUnitOperationError, operator.add, a1, a2
yield assert_raises, YTUnitOperationError, operator.iadd, a1, a2
+ # adding with zero is allowed irrespective of the units
+ zeros = np.zeros(3)
+ zeros_yta_dimless = YTArray(zeros, 'dimensionless')
+ zeros_yta_length = YTArray(zeros, 'm')
+ zeros_yta_mass = YTArray(zeros, 'kg')
+ operands = [0, YTQuantity(0), YTQuantity(0, 'kg'), zeros, zeros_yta_dimless,
+ zeros_yta_length, zeros_yta_mass]
+
+ for op in [operator.add, np.add]:
+ for operand in operands:
+ yield operate_and_compare, a1, operand, op, a1
+ yield operate_and_compare, operand, a1, op, a1
+ yield operate_and_compare, 4*m, operand, op, 4*m
+ yield operate_and_compare, operand, 4*m, op, 4*m
def test_subtraction():
"""
@@ -173,6 +187,20 @@
yield assert_raises, YTUnitOperationError, operator.sub, a1, a2
yield assert_raises, YTUnitOperationError, operator.isub, a1, a2
+ # subtracting with zero is allowed irrespective of the units
+ zeros = np.zeros(3)
+ zeros_yta_dimless = YTArray(zeros, 'dimensionless')
+ zeros_yta_length = YTArray(zeros, 'm')
+ zeros_yta_mass = YTArray(zeros, 'kg')
+ operands = [0, YTQuantity(0), YTQuantity(0, 'kg'), zeros, zeros_yta_dimless,
+ zeros_yta_length, zeros_yta_mass]
+
+ for op in [operator.sub, np.subtract]:
+ for operand in operands:
+ yield operate_and_compare, a1, operand, op, a1
+ yield operate_and_compare, operand, a1, op, -a1
+ yield operate_and_compare, 4*m, operand, op, 4*m
+ yield operate_and_compare, operand, 4*m, op, -4*m
def test_multiplication():
"""
@@ -1163,3 +1191,10 @@
degree_values = np.random.random(10)*degree
radian_values = degree_values.in_units('radian')
assert_array_equal(ufunc(degree_values), ufunc(radian_values))
+
+def test_builtin_sum():
+ from yt.units import km
+
+ arr = [1, 2, 3]*km
+ assert_equal(sum(arr), 6*km)
+
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/units/unit_object.py
--- a/yt/units/unit_object.py
+++ b/yt/units/unit_object.py
@@ -366,8 +366,12 @@
""" Test unit inequality. """
if not isinstance(u, Unit):
return True
- return \
- (self.base_value != u.base_value or self.dimensions != u.dimensions)
+ if self.base_value != u.base_value:
+ return True
+ # use 'is' comparison dimensions to avoid expensive sympy operation
+ if self.dimensions is u.dimensions:
+ return False
+ return self.dimensions != u.dimensions
def copy(self):
return copy.deepcopy(self)
@@ -389,6 +393,9 @@
def same_dimensions_as(self, other_unit):
""" Test if dimensions are the same. """
+ # test first for 'is' equality to avoid expensive sympy operation
+ if self.dimensions is other_unit.dimensions:
+ return True
return (self.dimensions / other_unit.dimensions) == sympy_one
@property
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -148,12 +148,18 @@
# attribute.
if isinstance(ret, YTArray):
if not inp.units.same_dimensions_as(ret.units):
+ # handle special case of adding or subtracting with zero or
+ # array filled with zero
+ if not np.any(other_object):
+ return ret.view(np.ndarray)
+ elif not np.any(this_object):
+ return ret
raise YTUnitOperationError(op_string, inp.units, ret.units)
ret = ret.in_units(inp.units)
- # If the other object is not a YTArray, the only valid case is adding
- # dimensionless things.
else:
- if not inp.units.is_dimensionless:
+ # If the other object is not a YTArray, then one of the arrays must be
+ # dimensionless or filled with zeros
+ if not inp.units.is_dimensionless and np.any(ret):
raise YTUnitOperationError(op_string, inp.units, dimensionless)
return ret
@@ -168,6 +174,31 @@
return other
+ at lru_cache(maxsize=128, typed=False)
+def _unit_repr_check_same(my_units, other_units):
+ """
+ Takes a Unit object, or string of known unit symbol, and check that it
+ is compatible with this quantity. Returns Unit object.
+
+ """
+ # let Unit() handle units arg if it's not already a Unit obj.
+ if not isinstance(other_units, Unit):
+ other_units = Unit(other_units, registry=my_units.registry)
+
+ equiv_dims = em_dimensions.get(my_units.dimensions, None)
+ if equiv_dims == other_units.dimensions:
+ if current_mks in equiv_dims.free_symbols:
+ base = "SI"
+ else:
+ base = "CGS"
+ raise YTEquivalentDimsError(my_units, other_units, base)
+
+ if not my_units.same_dimensions_as(other_units):
+ raise YTUnitConversionError(
+ my_units, my_units.dimensions, other_units, other_units.dimensions)
+
+ return other_units
+
unary_operators = (
negative, absolute, rint, ones_like, sign, conj, exp, exp2, log, log2,
log10, expm1, log1p, sqrt, square, reciprocal, sin, cos, tan, arcsin,
@@ -424,30 +455,6 @@
# Start unit conversion methods
#
- def _unit_repr_check_same(self, units):
- """
- Takes a Unit object, or string of known unit symbol, and check that it
- is compatible with this quantity. Returns Unit object.
-
- """
- # let Unit() handle units arg if it's not already a Unit obj.
- if not isinstance(units, Unit):
- units = Unit(units, registry=self.units.registry)
-
- equiv_dims = em_dimensions.get(self.units.dimensions,None)
- if equiv_dims == units.dimensions:
- if current_mks in equiv_dims.free_symbols:
- base = "SI"
- else:
- base = "CGS"
- raise YTEquivalentDimsError(self.units, units, base)
-
- if not self.units.same_dimensions_as(units):
- raise YTUnitConversionError(
- self.units, self.units.dimensions, units, units.dimensions)
-
- return units
-
def convert_to_units(self, units):
"""
Convert the array and units to the given units.
@@ -458,7 +465,7 @@
The units you want to convert to.
"""
- new_units = self._unit_repr_check_same(units)
+ new_units = _unit_repr_check_same(self.units, units)
(conversion_factor, offset) = self.units.get_conversion_factor(new_units)
self.units = new_units
@@ -517,11 +524,10 @@
YTArray
"""
- new_units = self._unit_repr_check_same(units)
+ new_units = _unit_repr_check_same(self.units, units)
(conversion_factor, offset) = self.units.get_conversion_factor(new_units)
- new_array = self * conversion_factor
- new_array.units = new_units
+ new_array = type(self)(self.ndview * conversion_factor, new_units)
if offset:
np.subtract(new_array, offset*new_array.uq, new_array)
@@ -1202,13 +1208,21 @@
unit2 = 1.0
unit_operator = self._ufunc_registry[context[0]]
if unit_operator in (preserve_units, comparison_unit, arctan2_unit):
- # Allow comparisons with dimensionless quantities.
- if (unit1 != unit2 and
- not unit2.is_dimensionless and not unit1.is_dimensionless):
- if not unit1.same_dimensions_as(unit2):
- raise YTUnitOperationError(context[0], unit1, unit2)
- else:
- raise YTUfuncUnitError(context[0], unit1, unit2)
+ # Allow comparisons, addition, and subtraction with
+ # dimensionless quantities or arrays filled with zeros.
+ u1d = unit1.is_dimensionless
+ u2d = unit2.is_dimensionless
+ if unit1 != unit2:
+ any_nonzero = [np.any(oper1), np.any(oper2)]
+ if any_nonzero[0] is np.bool_(False):
+ unit1 = unit2
+ elif any_nonzero[1] is np.bool_(False):
+ unit2 = unit1
+ elif not any([u1d, u2d]):
+ if not unit1.same_dimensions_as(unit2):
+ raise YTUnitOperationError(context[0], unit1, unit2)
+ else:
+ raise YTUfuncUnitError(context[0], unit1, unit2)
unit = unit_operator(unit1, unit2)
if unit_operator in (multiply_units, divide_units):
if unit.is_dimensionless and unit.base_value != 1.0:
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/utilities/command_line.py
--- a/yt/utilities/command_line.py
+++ b/yt/utilities/command_line.py
@@ -692,7 +692,10 @@
"\nHi there! Welcome to yt.\n\nWe've loaded your dataset as 'ds'. Enjoy!"
)
else:
- from IPython.config.loader import Config
+ try:
+ from traitlets.config.loader import Config
+ except ImportError:
+ from IPython.config.loader import Config
import sys
cfg = Config()
# prepend sys.path with current working directory
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/utilities/lib/amr_kdtools.pyx
--- a/yt/utilities/lib/amr_kdtools.pyx
+++ b/yt/utilities/lib/amr_kdtools.pyx
@@ -438,7 +438,6 @@
# If best_dim is -1, then we have found a place where there are no choices.
# Exit out and set the node to None.
if best_dim == -1:
- print 'Failed to split grid.'
return -1
@@ -509,6 +508,11 @@
best_dim = dim
my_max = n_unique
my_split = (n_unique-1)/2
+ if best_dim == -1:
+ for i in range(3):
+ free(uniquedims[i])
+ free(uniquedims)
+ return -1, 0, 0, 0
# I recognize how lame this is.
cdef np.ndarray[np.float64_t, ndim=1] tarr = np.empty(my_max, dtype='float64')
for i in range(my_max):
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/utilities/lib/mesh_construction.h
--- a/yt/utilities/lib/mesh_construction.h
+++ b/yt/utilities/lib/mesh_construction.h
@@ -37,7 +37,7 @@
{-1, -1, -1}
};
-// Triangule wedges
+// Triangulate wedges
int triangulate_wedge[MAX_NUM_TRI][3] = {
{0, 1, 2},
{0, 3, 1},
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/utilities/lib/primitives.pxd
--- a/yt/utilities/lib/primitives.pxd
+++ b/yt/utilities/lib/primitives.pxd
@@ -32,17 +32,17 @@
cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil
-cdef inline np.int64_t ray_triangle_intersect(const void* primitives,
- const np.int64_t item,
- Ray* ray) nogil
+cdef np.int64_t ray_triangle_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
-cdef inline void triangle_centroid(const void *primitives,
- const np.int64_t item,
- np.float64_t[3] centroid) nogil
+cdef void triangle_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
-cdef inline void triangle_bbox(const void *primitives,
- const np.int64_t item,
- BBox* bbox) nogil
+cdef void triangle_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
cdef struct Patch:
np.float64_t[8][3] v # 8 vertices per patch
@@ -67,14 +67,14 @@
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil
-cdef inline np.int64_t ray_patch_intersect(const void* primitives,
- const np.int64_t item,
- Ray* ray) nogil
+cdef np.int64_t ray_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
-cdef inline void patch_centroid(const void *primitives,
- const np.int64_t item,
- np.float64_t[3] centroid) nogil
+cdef void patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
-cdef inline void patch_bbox(const void *primitives,
- const np.int64_t item,
- BBox* bbox) nogil
+cdef void patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -35,9 +35,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline np.int64_t ray_triangle_intersect(const void* primitives,
- const np.int64_t item,
- Ray* ray) nogil:
+cdef np.int64_t ray_triangle_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
# https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
cdef Triangle tri = (<Triangle*> primitives)[item]
@@ -84,9 +84,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void triangle_centroid(const void *primitives,
- const np.int64_t item,
- np.float64_t[3] centroid) nogil:
+cdef void triangle_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
cdef Triangle tri = (<Triangle*> primitives)[item]
cdef np.int64_t i
@@ -97,9 +97,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void triangle_bbox(const void *primitives,
- const np.int64_t item,
- BBox* bbox) nogil:
+cdef void triangle_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
cdef Triangle tri = (<Triangle*> primitives)[item]
cdef np.int64_t i
@@ -240,9 +240,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline np.int64_t ray_patch_intersect(const void* primitives,
- const np.int64_t item,
- Ray* ray) nogil:
+cdef np.int64_t ray_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
cdef Patch patch = (<Patch*> primitives)[item]
@@ -264,9 +264,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void patch_centroid(const void *primitives,
- const np.int64_t item,
- np.float64_t[3] centroid) nogil:
+cdef void patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
cdef np.int64_t i, j
cdef Patch patch = (<Patch*> primitives)[item]
@@ -285,9 +285,9 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void patch_bbox(const void *primitives,
- const np.int64_t item,
- BBox* bbox) nogil:
+cdef void patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
cdef np.int64_t i, j
cdef Patch patch = (<Patch*> primitives)[item]
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/visualization/volume_rendering/interactive_vr.py
--- a/yt/visualization/volume_rendering/interactive_vr.py
+++ b/yt/visualization/volume_rendering/interactive_vr.py
@@ -18,6 +18,7 @@
from collections import OrderedDict
import matplotlib.cm as cm
import numpy as np
+import ctypes
from yt.utilities.math_utils import \
get_translate_matrix, \
@@ -29,7 +30,6 @@
rotation_matrix_to_quaternion
from .shader_objects import known_shaders, ShaderProgram
-
bbox_vertices = np.array(
[[ 0., 0., 0., 1.],
[ 0., 0., 1., 1.],
@@ -77,6 +77,27 @@
+1.0, +1.0, 0.0], dtype=np.float32
)
+triangulate_hex = np.array([
+ [0, 2, 1], [0, 3, 2],
+ [4, 5, 6], [4, 6, 7],
+ [0, 1, 5], [0, 5, 4],
+ [1, 2, 6], [1, 6, 5],
+ [0, 7, 3], [0, 4, 7],
+ [3, 6, 2], [3, 7, 6]]
+)
+
+triangulate_tetra = np.array([
+ [0, 1, 3], [2, 3, 1],
+ [0, 3, 2], [0, 2, 1]]
+)
+
+triangulate_wedge = np.array([
+ [3, 0, 1], [4, 3, 1],
+ [2, 5, 4], [2, 4, 1],
+ [0, 3, 2], [2, 3, 5],
+ [3, 4, 5], [0, 2, 1]]
+)
+
class IDVCamera(object):
'''Camera object used in the Interactive Data Visualization
@@ -532,8 +553,180 @@
GL.GL_RED, GL.GL_FLOAT, n_data.T)
GL.glGenerateMipmap(GL.GL_TEXTURE_3D)
+class ColorBarSceneComponent(SceneComponent):
+ '''
-class SceneGraph(SceneComponent):
+ A class for scene components that apply colorbars using a 1D texture.
+
+ '''
+
+ def __init__(self):
+ super(ColorBarSceneComponent, self).__init__()
+ self.camera = None
+ self.cmap_texture = None
+
+ def set_camera(self, camera):
+ pass
+
+ def update_minmax(self):
+ pass
+
+ def setup_cmap_tex(self):
+ '''Creates 1D texture that will hold colormap in framebuffer'''
+ self.cmap_texture = GL.glGenTextures(1) # create target texture
+ GL.glBindTexture(GL.GL_TEXTURE_1D, self.cmap_texture)
+ GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT, 1)
+ GL.glTexParameterf(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_WRAP_S, GL.GL_CLAMP_TO_EDGE)
+ GL.glTexParameteri(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_MAG_FILTER, GL.GL_LINEAR)
+ GL.glTexParameteri(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_MIN_FILTER, GL.GL_LINEAR)
+ GL.glTexImage1D(GL.GL_TEXTURE_1D, 0, GL.GL_RGBA, 256,
+ 0, GL.GL_RGBA, GL.GL_FLOAT, self.camera.cmap)
+ GL.glBindTexture(GL.GL_TEXTURE_1D, 0)
+
+ def update_cmap_tex(self):
+ '''Updates 1D texture with colormap that's used in framebuffer'''
+ if self.camera is None or not self.camera.cmap_new:
+ return
+
+ if self.cmap_texture is None:
+ self.setup_cmap_tex()
+
+ GL.glBindTexture(GL.GL_TEXTURE_1D, self.cmap_texture)
+ GL.glTexSubImage1D(GL.GL_TEXTURE_1D, 0, 0, 256,
+ GL.GL_RGBA, GL.GL_FLOAT, self.camera.cmap)
+ self.camera.cmap_new = False
+
+class MeshSceneComponent(ColorBarSceneComponent):
+ '''
+
+ A scene component for representing unstructured mesh data.
+
+ '''
+
+ def __init__(self, data_source, field):
+ super(MeshSceneComponent, self).__init__()
+ self.set_shader("mesh.v")
+ self.set_shader("mesh.f")
+
+ self.data_source = None
+ self.redraw = True
+
+ GL.glEnable(GL.GL_DEPTH_TEST)
+ GL.glDepthFunc(GL.GL_LESS)
+ GL.glEnable(GL.GL_CULL_FACE)
+ GL.glCullFace(GL.GL_BACK)
+
+ vertices, data, indices = self.get_mesh_data(data_source, field)
+
+ self._initialize_vertex_array("mesh_info")
+ GL.glBindVertexArray(self.vert_arrays["mesh_info"])
+
+ self.add_vert_attrib("vertex_buffer", vertices, vertices.size)
+ self.add_vert_attrib("data_buffer", data, data.size)
+
+ self.vert_attrib["element_buffer"] = (GL.glGenBuffers(1), indices.size)
+ GL.glBindBuffer(GL.GL_ELEMENT_ARRAY_BUFFER, self.vert_attrib["element_buffer"][0])
+ GL.glBufferData(GL.GL_ELEMENT_ARRAY_BUFFER, indices.nbytes, indices, GL.GL_STATIC_DRAW)
+
+ self.transform_matrix = GL.glGetUniformLocation(self.program.program,
+ "model_to_clip")
+
+ self.cmin = data.min()
+ self.cmax = data.max()
+
+ def set_camera(self, camera):
+ r""" Sets the camera orientation for the entire scene.
+
+ Parameters
+ ----------
+ camera : Camera
+
+ """
+ self.camera = camera
+ self.camera.cmap_min = float(self.cmin)
+ self.camera.cmap_max = float(self.cmax)
+ self.redraw = True
+
+ def get_mesh_data(self, data_source, field):
+ """
+
+ This reads the mesh data into a form that can be fed in to OpenGL.
+
+ """
+
+ # get mesh information
+ ftype, fname = field
+ mesh_id = int(ftype[-1])
+ mesh = data_source.ds.index.meshes[mesh_id-1]
+ offset = mesh._index_offset
+ vertices = mesh.connectivity_coords
+ indices = mesh.connectivity_indices - offset
+
+ # get vertex data
+ data = data_source[field]
+ vertex_data = np.zeros(vertices.shape[0], dtype=data.dtype)
+ vertex_data[indices.flatten()] = data.flatten()
+
+ if indices.shape[1] == 8:
+ tri_array = triangulate_hex
+ elif indices.shape[1] == 4:
+ tri_array = triangulate_tetra
+ elif indices.shape[1] == 6:
+ tri_array = triangulate_wedge
+ else:
+ raise NotImplementedError
+
+ tri_indices = []
+ for elem in indices:
+ for tri in tri_array:
+ tri_indices.append(elem[tri])
+ tri_indices = np.array(tri_indices)
+
+ v = vertices.astype(np.float32).flatten()
+ d = vertex_data.astype(np.float32).flatten()
+ i = tri_indices.astype(np.uint32).flatten()
+
+ return v, d, i
+
+ def run_program(self):
+ """ Renders one frame of the scene. """
+ with self.program.enable():
+
+ # Handle colormap
+ self.update_cmap_tex()
+
+ GL.glClear(GL.GL_COLOR_BUFFER_BIT | GL.GL_DEPTH_BUFFER_BIT)
+ projection_matrix = self.camera.projection_matrix
+ view_matrix = self.camera.view_matrix
+ model_to_clip = np.dot(projection_matrix, view_matrix)
+ GL.glUniformMatrix4fv(self.transform_matrix, 1, True, model_to_clip)
+
+ GL.glActiveTexture(GL.GL_TEXTURE1)
+ GL.glBindTexture(GL.GL_TEXTURE_1D, self.cmap_texture)
+
+ self.program._set_uniform("cmap", 0)
+ self.program._set_uniform("cmap_min", self.camera.cmap_min)
+ self.program._set_uniform("cmap_max", self.camera.cmap_max)
+
+ GL.glEnableVertexAttribArray(0)
+ GL.glBindBuffer(GL.GL_ARRAY_BUFFER, self.vert_attrib["vertex_buffer"][0])
+ GL.glVertexAttribPointer(0, 3, GL.GL_FLOAT, False, 0, ctypes.c_void_p(0))
+
+ GL.glEnableVertexAttribArray(1)
+ GL.glBindBuffer(GL.GL_ARRAY_BUFFER, self.vert_attrib["data_buffer"][0])
+ GL.glVertexAttribPointer(1, 1, GL.GL_FLOAT, False, 0, ctypes.c_void_p(0))
+
+ GL.glBindBuffer(GL.GL_ELEMENT_ARRAY_BUFFER, self.vert_attrib["element_buffer"][0])
+ GL.glDrawElements(GL.GL_TRIANGLES, self.vert_attrib["element_buffer"][1],
+ GL.GL_UNSIGNED_INT, ctypes.c_void_p(0))
+
+ GL.glDisableVertexAttribArray(0)
+ GL.glDisableVertexAttribArray(1)
+
+ render = run_program
+
+
+class SceneGraph(ColorBarSceneComponent):
"""A basic OpenGL render for IDV.
The SceneGraph class is the primary driver behind creating a IDV rendering.
@@ -554,8 +747,6 @@
self.collections = []
self.fbo = None
self.fb_texture = None
- self.cmap_texture = None
- self.camera = None
self.shader_program = None
self.fb_shader_program = None
self.min_val, self.max_val = 1e60, -1e60
@@ -587,36 +778,8 @@
self.setup_fb(self.width, self.height)
- def setup_cmap_tex(self):
- '''Creates 1D texture that will hold colormap in framebuffer'''
- self.cmap_texture = GL.glGenTextures(1) # create target texture
- GL.glBindTexture(GL.GL_TEXTURE_1D, self.cmap_texture)
- GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT, 1)
- GL.glTexParameterf(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_WRAP_S, GL.GL_CLAMP_TO_EDGE)
- GL.glTexParameteri(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_MAG_FILTER, GL.GL_LINEAR)
- GL.glTexParameteri(GL.GL_TEXTURE_1D, GL.GL_TEXTURE_MIN_FILTER, GL.GL_LINEAR)
- GL.glTexImage1D(GL.GL_TEXTURE_1D, 0, GL.GL_RGBA, 256,
- 0, GL.GL_RGBA, GL.GL_FLOAT, self.camera.cmap)
- GL.glBindTexture(GL.GL_TEXTURE_1D, 0)
-
-
- def update_cmap_tex(self):
- '''Updates 1D texture with colormap that's used in framebuffer'''
- if self.camera is None or not self.camera.cmap_new:
- return
-
- if self.cmap_texture is None:
- self.setup_cmap_tex()
-
- GL.glBindTexture(GL.GL_TEXTURE_1D, self.cmap_texture)
- GL.glTexSubImage1D(GL.GL_TEXTURE_1D, 0, 0, 256,
- GL.GL_RGBA, GL.GL_FLOAT, self.camera.cmap)
- GL.glBindTexture(GL.GL_TEXTURE_1D, 0)
- self.camera.cmap_new = False
-
-
def setup_fb(self, width, height):
- '''Setups FrameBuffer that will be used as container
+ '''Sets up FrameBuffer that will be used as container
for 1 pass of rendering'''
# Clean up old FB and Texture
if self.fb_texture is not None and \
@@ -670,7 +833,6 @@
status = GL.glCheckFramebufferStatus(GL.GL_FRAMEBUFFER)
assert status == GL.GL_FRAMEBUFFER_COMPLETE, status
-
def add_collection(self, collection):
r"""Adds a block collection to the scene. Collections must not overlap.
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/visualization/volume_rendering/interactive_vr_helpers.py
--- a/yt/visualization/volume_rendering/interactive_vr_helpers.py
+++ b/yt/visualization/volume_rendering/interactive_vr_helpers.py
@@ -58,7 +58,8 @@
raise ImportError("This functionality requires the cyglfw3 and PyOpenGL "
"packages to be installed.")
- from .interactive_vr import SceneGraph, BlockCollection, TrackballCamera
+ from .interactive_vr import SceneGraph, BlockCollection, TrackballCamera, \
+ MeshSceneComponent
from .interactive_loop import RenderingContext
if isinstance(data_source, Dataset):
@@ -78,16 +79,23 @@
if cam_focus is None:
cam_focus = dobj.ds.domain_center
+ rc = RenderingContext(*window_size)
+
+ if hasattr(dobj.ds.index, "meshes"):
+ # unstructured mesh datasets tend to have tight
+ # domain boundaries, do some extra padding here.
+ cam_position = 3.0*dobj.ds.domain_right_edge
+ scene = MeshSceneComponent(dobj, field)
+ else:
+ scene = SceneGraph()
+ collection = BlockCollection()
+ collection.add_data(dobj, field)
+ scene.add_collection(collection)
+
aspect_ratio = window_size[1] / window_size[0]
far_plane = np.linalg.norm(cam_focus - cam_position) * 2.0
near_plane = 0.01 * far_plane
- rc = RenderingContext(*window_size)
- scene = SceneGraph()
- collection = BlockCollection()
- collection.add_data(dobj, field)
- scene.add_collection(collection)
-
c = TrackballCamera(position=cam_position, focus=cam_focus, near_plane=near_plane,
far_plane=far_plane, aspect_ratio=aspect_ratio)
rc.start_loop(scene, c)
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/visualization/volume_rendering/shader_objects.py
--- a/yt/visualization/volume_rendering/shader_objects.py
+++ b/yt/visualization/volume_rendering/shader_objects.py
@@ -29,7 +29,7 @@
class ShaderProgram(object):
'''
Wrapper class that compiles and links vertex and fragment shaders
- into shader program.
+ into a shader program.
Parameters
----------
@@ -269,3 +269,13 @@
'''A second pass vertex shader that performs no operations on vertices'''
_source = "passthrough.vertexshader"
_shader_name = "passthrough.v"
+
+class MeshVertexShader(VertexShader):
+ '''A vertex shader used for unstructured mesh rendering.'''
+ _source = "mesh.vertexshader"
+ _shader_name = "mesh.v"
+
+class MeshFragmentShader(FragmentShader):
+ '''A vertex shader used for unstructured mesh rendering.'''
+ _source = "mesh.fragmentshader"
+ _shader_name = "mesh.f"
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/visualization/volume_rendering/shaders/mesh.fragmentshader
--- /dev/null
+++ b/yt/visualization/volume_rendering/shaders/mesh.fragmentshader
@@ -0,0 +1,17 @@
+#version 330 core
+
+in float fragmentData;
+out vec4 color;
+
+uniform sampler1D cmap;
+uniform float cmap_min;
+uniform float cmap_max;
+
+void main()
+{
+ float data = fragmentData;
+ float cm = cmap_min;
+ float cp = cmap_max;
+
+ color = texture(cmap, (data - cm) / (cp - cm));
+}
diff -r eb9e446f13191994fe3b92b6bfea64315504ebbe -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 yt/visualization/volume_rendering/shaders/mesh.vertexshader
--- /dev/null
+++ b/yt/visualization/volume_rendering/shaders/mesh.vertexshader
@@ -0,0 +1,11 @@
+#version 330 core
+
+layout(location = 0) in vec3 vertexPosition_modelspace;
+layout(location = 1) in float vertexData;
+out float fragmentData;
+uniform mat4 model_to_clip;
+void main()
+{
+ gl_Position = model_to_clip * vec4(vertexPosition_modelspace, 1);
+ fragmentData = vertexData;
+}
https://bitbucket.org/yt_analysis/yt/commits/d39bece014c4/
Changeset: d39bece014c4
Branch: yt
User: atmyers
Date: 2016-05-23 19:47:35+00:00
Summary: Move mesh triangulation code to the Cython utilities. Faster and avoids duplicating the various triangulation arrays.
Affected #: 1 file
diff -r 3bd9cbb7594b8994fee8fd1d7a63a8d7b6732230 -r d39bece014c4649162b7f378cb1d27be8c35842b yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -56,6 +56,38 @@
int hex20_faces[6][8]
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
+ cdef np.int64_t num_elem = indices.shape[0]
+ cdef np.int64_t num_verts = indices.shape[1]
+ cdef int[MAX_NUM_TRI][3] tri_array
+ cdef np.int64_t tri_per_elem
+
+ if (num_verts == 8 or num_verts == 20 or num_verts == 27):
+ tri_per_elem = HEX_NT
+ tri_array = triangulate_hex
+ elif num_verts == 4:
+ tri_per_elem = TETRA_NT
+ tri_array = triangulate_tetra
+ elif num_verts == 6:
+ tri_per_elem = WEDGE_NT
+ tri_array = triangulate_wedge
+ else:
+ raise YTElementTypeNotRecognized(3, num_verts)
+
+ cdef np.ndarray[np.int64_t, ndim=2] tri_indices
+ tri_indices = np.empty((tri_per_elem * num_elem, 3), dtype="int64")
+
+ cdef np.int64_t i, j
+ for i in range(num_elem):
+ for j in range(tri_per_elem):
+ for k in range(3):
+ tri_indices[i*tri_per_elem + j][k] = indices[i][tri_array[j][k]]
+ return tri_indices
+
+
cdef class LinearElementMesh:
r'''
https://bitbucket.org/yt_analysis/yt/commits/abc633ee1ce0/
Changeset: abc633ee1ce0
Branch: yt
User: atmyers
Date: 2016-05-23 19:48:07+00:00
Summary: Loop over all the meshes unless only one is requested.
Affected #: 1 file
diff -r d39bece014c4649162b7f378cb1d27be8c35842b -r abc633ee1ce03e85279811bc83063f57f26f72cf yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -35,31 +35,11 @@
periodic_ray
from yt.utilities.lib.line_integral_convolution import \
line_integral_convolution_2d
-from yt.utilities.exceptions import YTElementTypeNotRecognized
+from yt.geometry.unstructured_mesh_handler import UnstructuredIndex
+from yt.utilities.lib.mesh_construction import triangulate_indices
from . import _MPL
-triangulate_hex = np.array([
- [0, 2, 1], [0, 3, 2],
- [4, 5, 6], [4, 6, 7],
- [0, 1, 5], [0, 5, 4],
- [1, 2, 6], [1, 6, 5],
- [0, 7, 3], [0, 4, 7],
- [3, 6, 2], [3, 7, 6]]
-)
-
-triangulate_tetra = np.array([
- [0, 1, 3], [2, 3, 1],
- [0, 3, 2], [0, 2, 1]]
-)
-
-triangulate_wedge = np.array([
- [3, 0, 1], [4, 3, 1],
- [2, 5, 4], [2, 4, 1],
- [0, 3, 2], [2, 3, 5],
- [3, 4, 5], [0, 2, 1]]
-)
-
callback_registry = {}
class RegisteredCallback(type):
@@ -1626,7 +1606,9 @@
"""
annotate_mesh_lines()
- Adds the mesh lines to the plot.
+ Adds mesh lines to the plot. Only works for unstructured or
+ semi-structured mesh data. For structured grid data, see
+ GridBoundaryCallback or CellEdgesCallback.
Parameters
----------
@@ -1664,44 +1646,57 @@
return new_coords, new_connects
+ def _get_mesh(self, plot):
+
+ index = plot.ds.index
+ if not issubclass(type(index), UnstructuredIndex):
+ raise RuntimeError("Mesh line annotations only work for "
+ "unstructured or semi-structured mesh data.")
+
+ ftype, fname = plot.field
+
+ if ftype.startswith('connect'):
+ # select appropriate mesh from file
+ mesh_id = int(ftype[-1]) - 1
+ mesh = index.meshes[mesh_id]
+ coords = mesh.connectivity_coords
+ indices = mesh.connectivity_indices - mesh._index_offset
+ else:
+ # select all the meshes
+ m = index.meshes[0]
+ coords = m.connectivity_coords
+ indices = m.connectivity_indices - m._index_offset
+ for m in index.meshes[1:]:
+ next_indices = m.connectivity_indices - m._index_offset
+ indices = np.vstack((indices, next_indices))
+ return coords, indices
+
def __call__(self, plot):
index = plot.ds.index
+ if not issubclass(type(index), UnstructuredIndex):
+ raise RuntimeError("Mesh line annotations only work for "
+ "unstructured or semi-structured mesh data.")
ftype, fname = plot.field
- mesh_id = int(ftype[-1]) - 1
- mesh = index.meshes[mesh_id]
- coords = mesh.connectivity_coords
- indices = mesh.connectivity_indices - mesh._index_offset
+ for i, m in enumerate(index.meshes):
+ if ftype.startswith('connect') and int(ftype[-1]) - 1 != i:
+ continue
+ coords = m.connectivity_coords
+ indices = m.connectivity_indices - m._index_offset
+
+ num_verts = indices.shape[1]
+ num_dims = coords.shape[1]
- num_verts = indices.shape[1]
- num_dims = coords.shape[1]
+ if num_dims == 2 and num_verts == 3:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
+ elif num_dims == 2 and num_verts == 4:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
- if num_dims == 3 and (num_verts == 8 or
- num_verts == 20 or
- num_verts == 27):
- tri_array = triangulate_hex
- elif num_dims == 3 and num_verts == 4:
- tri_array = triangulate_tetra
- elif num_dims == 3 and num_verts == 6:
- tri_array = triangulate_wedge
- elif num_dims == 2 and num_verts == 3:
- tri_array = triangulate_wedge
- coords, indices = self.promote_2d_to_3d(coords, indices, plot)
- elif num_dims == 2 and num_verts == 4:
- tri_array = triangulate_hex
- coords, indices = self.promote_2d_to_3d(coords, indices, plot)
- else:
- raise YTElementTypeNotRecognized(num_dims, num_verts)
-
- tri_indices = []
- for elem in indices:
- for tri in tri_array:
- tri_indices.append(elem[tri])
- tri_indices = np.array(tri_indices)
- points = coords[tri_indices]
+ tri_indices = triangulate_indices(indices)
+ points = coords[tri_indices]
- tfc = TriangleFacetsCallback(points, plot_args=self.plot_args)
- return tfc(plot)
+ tfc = TriangleFacetsCallback(points, plot_args=self.plot_args)
+ tfc(plot)
class TriangleFacetsCallback(PlotCallback):
https://bitbucket.org/yt_analysis/yt/commits/e1c1755627ec/
Changeset: e1c1755627ec
Branch: yt
User: atmyers
Date: 2016-05-23 19:50:57+00:00
Summary: Use a regular RuntimeError here, as suggested by Nathan.
Affected #: 1 file
diff -r abc633ee1ce03e85279811bc83063f57f26f72cf -r e1c1755627ec1ddce93f2ea3ca6489016de0aea9 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -20,7 +20,6 @@
from yt.utilities.lib.fp_utils cimport fmin, fmax, i64min, i64max, imin, imax, fabs
from yt.utilities.exceptions import \
YTPixelizeError, \
- YTException, \
YTElementTypeNotRecognized
from libc.stdlib cimport malloc, free
from vec3_ops cimport dot, cross, subtract
@@ -593,10 +592,9 @@
# if we are in 2D land, the 1 cell thick dimension had better be 'z'
if ndim == 2:
- try:
- assert(buff_size[2] == 1)
- except AssertionError:
- raise YTException("2D datasets can only be sliced in the 'z' direction.")
+ if buff_size[2] != 1:
+ raise RuntimeError("Slices of 2D datasets must be "
+ "perpendicular to the 'z' direction.")
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
https://bitbucket.org/yt_analysis/yt/commits/fecd15177d5b/
Changeset: fecd15177d5b
Branch: yt
User: atmyers
Date: 2016-05-23 20:02:31+00:00
Summary: adding a test of the mesh lines callback to test_callbacks.py
Affected #: 1 file
diff -r e1c1755627ec1ddce93f2ea3ca6489016de0aea9 -r fecd15177d5b844519a04c450df9df8d8b3b9dc3 yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -21,7 +21,9 @@
from yt.config import \
ytcfg
from yt.testing import \
- fake_amr_ds
+ fake_amr_ds, \
+ fake_tetrahedral_ds, \
+ fake_hexahedral_ds
import yt.units as u
from .test_plotwindow import assert_fname
from yt.utilities.exceptions import \
@@ -368,6 +370,22 @@
color=(0.0, 1.0, 1.0))
p.save(prefix)
+def test_mesh_lines_callback():
+ with _cleanup_fname() as prefix:
+
+ ds = fake_hexahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+ ds = fake_tetrahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+
def test_line_integral_convolution_callback():
with _cleanup_fname() as prefix:
ds = fake_amr_ds(fields =
@@ -390,3 +408,4 @@
alpha=0.9, const_alpha=True)
p.save(prefix)
+
https://bitbucket.org/yt_analysis/yt/commits/6a6e2e7ebc74/
Changeset: 6a6e2e7ebc74
Branch: yt
User: atmyers
Date: 2016-05-23 20:04:20+00:00
Summary: remove now unused _get_mesh method.
Affected #: 1 file
diff -r fecd15177d5b844519a04c450df9df8d8b3b9dc3 -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -1646,31 +1646,6 @@
return new_coords, new_connects
- def _get_mesh(self, plot):
-
- index = plot.ds.index
- if not issubclass(type(index), UnstructuredIndex):
- raise RuntimeError("Mesh line annotations only work for "
- "unstructured or semi-structured mesh data.")
-
- ftype, fname = plot.field
-
- if ftype.startswith('connect'):
- # select appropriate mesh from file
- mesh_id = int(ftype[-1]) - 1
- mesh = index.meshes[mesh_id]
- coords = mesh.connectivity_coords
- indices = mesh.connectivity_indices - mesh._index_offset
- else:
- # select all the meshes
- m = index.meshes[0]
- coords = m.connectivity_coords
- indices = m.connectivity_indices - m._index_offset
- for m in index.meshes[1:]:
- next_indices = m.connectivity_indices - m._index_offset
- indices = np.vstack((indices, next_indices))
- return coords, indices
-
def __call__(self, plot):
index = plot.ds.index
https://bitbucket.org/yt_analysis/yt/commits/413338364456/
Changeset: 413338364456
Branch: yt
User: atmyers
Date: 2016-05-23 21:42:00+00:00
Summary: some refactoring to allow the use of triangulate_indices with having embree installed.
Affected #: 6 files
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 setup.py
--- a/setup.py
+++ b/setup.py
@@ -157,12 +157,6 @@
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),
- Extension("yt.utilities.lib.primitives",
- ["yt/utilities/lib/primitives.pyx"],
- libraries=std_libs,
- depends=["yt/utilities/lib/primitives.pxd",
- "yt/utilities/lib/vec3_ops.pxd",
- "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.bounding_volume_hierarchy",
["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
include_dirs=["yt/utilities/lib/"],
@@ -170,6 +164,7 @@
extra_link_args=omp_args,
libraries=std_libs,
depends=["yt/utilities/lib/element_mappings.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/vec3_ops.pxd",
"yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
@@ -196,6 +191,9 @@
"yt/utilities/lib/fixed_interpolator.pxd",
"yt/utilities/lib/fixed_interpolator.h",
]),
+ Extension("yt.utilities.lib.mesh_triangulation",
+ ["yt/utilities/lib/mesh_triangulation.pyx"],
+ depends=["yt/utilities/lib/mesh_triangulation.h"]),
Extension("yt.utilities.lib.pixelization_routines",
["yt/utilities/lib/pixelization_routines.pyx",
"yt/utilities/lib/pixelization_constants.c"],
@@ -203,6 +201,12 @@
libraries=std_libs, depends=["yt/utilities/lib/fp_utils.pxd",
"yt/utilities/lib/pixelization_constants.h",
"yt/utilities/lib/element_mappings.pxd"]),
+ Extension("yt.utilities.lib.primitives",
+ ["yt/utilities/lib/primitives.pyx"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/primitives.pxd",
+ "yt/utilities/lib/vec3_ops.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.origami",
["yt/utilities/lib/origami.pyx",
"yt/utilities/lib/origami_tags.c"],
@@ -282,6 +286,7 @@
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],
depends=["yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/mesh_intersection.pxd",
"yt/utlilites/lib/mesh_samplers.pxd",
"yt/utlilites/lib/mesh_traversal.pxd"]),
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -4,7 +4,7 @@
from yt.utilities.lib.element_mappings cimport ElementSampler
from yt.utilities.lib.primitives cimport BBox, Ray
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -40,7 +40,7 @@
patchBoundsFunc
from yt.utilities.exceptions import YTElementTypeNotRecognized
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
@@ -56,38 +56,6 @@
int hex20_faces[6][8]
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
- cdef np.int64_t num_elem = indices.shape[0]
- cdef np.int64_t num_verts = indices.shape[1]
- cdef int[MAX_NUM_TRI][3] tri_array
- cdef np.int64_t tri_per_elem
-
- if (num_verts == 8 or num_verts == 20 or num_verts == 27):
- tri_per_elem = HEX_NT
- tri_array = triangulate_hex
- elif num_verts == 4:
- tri_per_elem = TETRA_NT
- tri_array = triangulate_tetra
- elif num_verts == 6:
- tri_per_elem = WEDGE_NT
- tri_array = triangulate_wedge
- else:
- raise YTElementTypeNotRecognized(3, num_verts)
-
- cdef np.ndarray[np.int64_t, ndim=2] tri_indices
- tri_indices = np.empty((tri_per_elem * num_elem, 3), dtype="int64")
-
- cdef np.int64_t i, j
- for i in range(num_elem):
- for j in range(tri_per_elem):
- for k in range(3):
- tri_indices[i*tri_per_elem + j][k] = indices[i][tri_array[j][k]]
- return tri_indices
-
-
cdef class LinearElementMesh:
r'''
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 yt/utilities/lib/mesh_triangulation.h
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -0,0 +1,66 @@
+#define MAX_NUM_TRI 12
+#define HEX_NV 8
+#define HEX_NT 12
+#define TETRA_NV 4
+#define TETRA_NT 4
+#define WEDGE_NV 6
+#define WEDGE_NT 8
+
+// This array is used to triangulate the hexahedral mesh elements
+// Each element has six faces with two triangles each.
+// The vertex ordering convention is assumed to follow that used
+// here: http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
+// Note that this is the case for Exodus II data.
+int triangulate_hex[MAX_NUM_TRI][3] = {
+ {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
+ {4, 5, 6}, {4, 6, 7}, // Face is 4 5 6 7
+ {0, 1, 5}, {0, 5, 4}, // Face is 0 1 5 4
+ {1, 2, 6}, {1, 6, 5}, // Face is 1 2 6 5
+ {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
+ {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
+};
+
+// Similarly, this is used to triangulate the tetrahedral cells
+int triangulate_tetra[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 1, 3},
+ {0, 2, 3},
+ {1, 2, 3},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+// Triangulate wedges
+int triangulate_wedge[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 3, 1},
+ {1, 3, 4},
+ {0, 2, 3},
+ {2, 5, 3},
+ {1, 4, 2},
+ {2, 4, 5},
+ {3, 5, 4},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+
+// This is used to select faces from a 20-sided hex element
+int hex20_faces[6][8] = {
+ {0, 1, 5, 4, 12, 8, 13, 16},
+ {1, 2, 6, 5, 13, 9, 14, 17},
+ {3, 2, 6, 7, 15, 10, 14, 18},
+ {0, 3, 7, 4, 12, 11, 15, 19},
+ {4, 5, 6, 7, 19, 16, 17, 18},
+ {0, 1, 2, 3, 11, 8, 9, 10}
+};
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 yt/utilities/lib/mesh_triangulation.pyx
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -0,0 +1,52 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+
+from yt.utilities.exceptions import YTElementTypeNotRecognized
+
+cdef extern from "mesh_triangulation.h":
+ enum:
+ MAX_NUM_TRI
+
+ int HEX_NV
+ int HEX_NT
+ int TETRA_NV
+ int TETRA_NT
+ int WEDGE_NV
+ int WEDGE_NT
+ int triangulate_hex[MAX_NUM_TRI][3]
+ int triangulate_tetra[MAX_NUM_TRI][3]
+ int triangulate_wedge[MAX_NUM_TRI][3]
+ int hex20_faces[6][8]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
+ cdef np.int64_t num_elem = indices.shape[0]
+ cdef np.int64_t num_verts = indices.shape[1]
+ cdef int[MAX_NUM_TRI][3] tri_array
+ cdef np.int64_t tri_per_elem
+
+ if (num_verts == 8 or num_verts == 20 or num_verts == 27):
+ tri_per_elem = HEX_NT
+ tri_array = triangulate_hex
+ elif num_verts == 4:
+ tri_per_elem = TETRA_NT
+ tri_array = triangulate_tetra
+ elif num_verts == 6:
+ tri_per_elem = WEDGE_NT
+ tri_array = triangulate_wedge
+ else:
+ raise YTElementTypeNotRecognized(3, num_verts)
+
+ cdef np.ndarray[np.int64_t, ndim=2] tri_indices
+ tri_indices = np.empty((tri_per_elem * num_elem, 3), dtype="int64")
+
+ cdef np.int64_t i, j
+ for i in range(num_elem):
+ for j in range(tri_per_elem):
+ for k in range(3):
+ tri_indices[i*tri_per_elem + j][k] = indices[i][tri_array[j][k]]
+ return tri_indices
diff -r 6a6e2e7ebc74bc238c689890b04f5da34a9d0728 -r 413338364456ba3ce6603427cd30a2dc46cebea8 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -36,7 +36,7 @@
from yt.utilities.lib.line_integral_convolution import \
line_integral_convolution_2d
from yt.geometry.unstructured_mesh_handler import UnstructuredIndex
-from yt.utilities.lib.mesh_construction import triangulate_indices
+from yt.utilities.lib.mesh_triangulation import triangulate_indices
from . import _MPL
https://bitbucket.org/yt_analysis/yt/commits/2b7d57da527f/
Changeset: 2b7d57da527f
Branch: yt
User: atmyers
Date: 2016-05-23 21:53:06+00:00
Summary: removing old file.
Affected #: 1 file
diff -r 413338364456ba3ce6603427cd30a2dc46cebea8 -r 2b7d57da527fc2cfd3d02248ebe1d0305cc12a99 yt/utilities/lib/mesh_construction.h
--- a/yt/utilities/lib/mesh_construction.h
+++ /dev/null
@@ -1,66 +0,0 @@
-#define MAX_NUM_TRI 12
-#define HEX_NV 8
-#define HEX_NT 12
-#define TETRA_NV 4
-#define TETRA_NT 4
-#define WEDGE_NV 6
-#define WEDGE_NT 8
-
-// This array is used to triangulate the hexahedral mesh elements
-// Each element has six faces with two triangles each.
-// The vertex ordering convention is assumed to follow that used
-// here: http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-// Note that this is the case for Exodus II data.
-int triangulate_hex[MAX_NUM_TRI][3] = {
- {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
- {4, 5, 6}, {4, 6, 7}, // Face is 4 5 6 7
- {0, 1, 5}, {0, 5, 4}, // Face is 0 1 5 4
- {1, 2, 6}, {1, 6, 5}, // Face is 1 2 6 5
- {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
- {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
-};
-
-// Similarly, this is used to triangulate the tetrahedral cells
-int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 1, 3},
- {0, 2, 3},
- {1, 2, 3},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-// Triangulate wedges
-int triangulate_wedge[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 3, 1},
- {1, 3, 4},
- {0, 2, 3},
- {2, 5, 3},
- {1, 4, 2},
- {2, 4, 5},
- {3, 5, 4},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-
-// This is used to select faces from a 20-sided hex element
-int hex20_faces[6][8] = {
- {0, 1, 5, 4, 12, 8, 13, 16},
- {1, 2, 6, 5, 13, 9, 14, 17},
- {3, 2, 6, 7, 15, 10, 14, 18},
- {0, 3, 7, 4, 12, 11, 15, 19},
- {4, 5, 6, 7, 19, 16, 17, 18},
- {0, 1, 2, 3, 11, 8, 9, 10}
-};
https://bitbucket.org/yt_analysis/yt/commits/db7542b48315/
Changeset: db7542b48315
Branch: yt
User: atmyers
Date: 2016-05-23 22:10:28+00:00
Summary: do not require plot.field to be a tuple.
Affected #: 1 file
diff -r 2b7d57da527fc2cfd3d02248ebe1d0305cc12a99 -r db7542b48315fe2362c1e84244fac712432732d6 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -1652,10 +1652,13 @@
if not issubclass(type(index), UnstructuredIndex):
raise RuntimeError("Mesh line annotations only work for "
"unstructured or semi-structured mesh data.")
- ftype, fname = plot.field
for i, m in enumerate(index.meshes):
- if ftype.startswith('connect') and int(ftype[-1]) - 1 != i:
- continue
+ try:
+ ftype, fname = plot.field
+ if ftype.startswith('connect') and int(ftype[-1]) - 1 != i:
+ continue
+ except ValueError:
+ pass
coords = m.connectivity_coords
indices = m.connectivity_indices - m._index_offset
https://bitbucket.org/yt_analysis/yt/commits/cdea0aeff19b/
Changeset: cdea0aeff19b
Branch: yt
User: atmyers
Date: 2016-05-24 02:29:48+00:00
Summary: faster triangulation routine.
Affected #: 1 file
diff -r db7542b48315fe2362c1e84244fac712432732d6 -r cdea0aeff19bae093e4db2d133a91732c5d2f440 yt/utilities/lib/mesh_triangulation.pyx
--- a/yt/utilities/lib/mesh_triangulation.pyx
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -25,28 +25,34 @@
@cython.cdivision(True)
def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
cdef np.int64_t num_elem = indices.shape[0]
- cdef np.int64_t num_verts = indices.shape[1]
+ cdef np.int64_t VPE = indices.shape[1] # num verts per element
+ cdef np.int64_t TPE # num triangles per element
cdef int[MAX_NUM_TRI][3] tri_array
- cdef np.int64_t tri_per_elem
-
- if (num_verts == 8 or num_verts == 20 or num_verts == 27):
- tri_per_elem = HEX_NT
+
+ if (VPE == 8 or VPE == 20 or VPE == 27):
+ TPE = HEX_NT
tri_array = triangulate_hex
- elif num_verts == 4:
- tri_per_elem = TETRA_NT
+ elif VPE == 4:
+ TPE = TETRA_NT
tri_array = triangulate_tetra
- elif num_verts == 6:
- tri_per_elem = WEDGE_NT
+ elif VPE == 6:
+ TPE = WEDGE_NT
tri_array = triangulate_wedge
else:
- raise YTElementTypeNotRecognized(3, num_verts)
+ raise YTElementTypeNotRecognized(3, VPE)
+ cdef np.int64_t num_tri = TPE * num_elem
+
cdef np.ndarray[np.int64_t, ndim=2] tri_indices
- tri_indices = np.empty((tri_per_elem * num_elem, 3), dtype="int64")
+ tri_indices = np.empty((num_tri, 3), dtype="int64")
+
+ cdef np.int64_t *tri_indices_ptr = <np.int64_t*> tri_indices.data
+ cdef np.int64_t *indices_ptr = <np.int64_t*> indices.data
- cdef np.int64_t i, j
+ cdef np.int64_t i, j, k, offset
for i in range(num_elem):
- for j in range(tri_per_elem):
+ for j in range(TPE):
for k in range(3):
- tri_indices[i*tri_per_elem + j][k] = indices[i][tri_array[j][k]]
+ offset = tri_array[j][k]
+ tri_indices_ptr[i*TPE*3 + 3*j + k] = indices_ptr[i*VPE + offset]
return tri_indices
https://bitbucket.org/yt_analysis/yt/commits/19aead8b2225/
Changeset: 19aead8b2225
Branch: yt
User: ngoldbaum
Date: 2016-05-24 03:30:21+00:00
Summary: Merged in atmyers/yt (pull request #2174)
Changing the way mesh line annotations are performed.
Affected #: 14 files
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -406,16 +406,14 @@
import yt
ds = yt.load('MOOSE_sample_data/out.e-s010')
sl = yt.SlicePlot(ds, 'z', ('connect1', 'diffused'))
- sl.annotate_mesh_lines(thresh=0.1)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
sl.zoom(0.75)
sl.save()
-This annotation is performed by marking the pixels where the mapped coordinate is close
-to the element boundary. What counts as 'close' (in the mapped coordinate system) is
-determined by the ``thresh`` parameter, which can be varied to make the lines thicker or
-thinner.
+The ``plot_args`` parameter is a dictionary of keyword arguments that will be passed
+to matplotlib. It can be used to control the mesh line color, thickness, etc...
-The above example all involve 8-node hexahedral mesh elements. Here is another example from
+The above examples all involve 8-node hexahedral mesh elements. Here is another example from
a dataset that uses 6-node wedge elements:
.. python-script::
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 setup.py
--- a/setup.py
+++ b/setup.py
@@ -157,12 +157,6 @@
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),
- Extension("yt.utilities.lib.primitives",
- ["yt/utilities/lib/primitives.pyx"],
- libraries=std_libs,
- depends=["yt/utilities/lib/primitives.pxd",
- "yt/utilities/lib/vec3_ops.pxd",
- "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.bounding_volume_hierarchy",
["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
include_dirs=["yt/utilities/lib/"],
@@ -170,6 +164,7 @@
extra_link_args=omp_args,
libraries=std_libs,
depends=["yt/utilities/lib/element_mappings.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/vec3_ops.pxd",
"yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
@@ -196,14 +191,22 @@
"yt/utilities/lib/fixed_interpolator.pxd",
"yt/utilities/lib/fixed_interpolator.h",
]),
+ Extension("yt.utilities.lib.mesh_triangulation",
+ ["yt/utilities/lib/mesh_triangulation.pyx"],
+ depends=["yt/utilities/lib/mesh_triangulation.h"]),
Extension("yt.utilities.lib.pixelization_routines",
["yt/utilities/lib/pixelization_routines.pyx",
"yt/utilities/lib/pixelization_constants.c"],
include_dirs=["yt/utilities/lib/"],
- language="c++",
libraries=std_libs, depends=["yt/utilities/lib/fp_utils.pxd",
"yt/utilities/lib/pixelization_constants.h",
"yt/utilities/lib/element_mappings.pxd"]),
+ Extension("yt.utilities.lib.primitives",
+ ["yt/utilities/lib/primitives.pyx"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/primitives.pxd",
+ "yt/utilities/lib/vec3_ops.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.origami",
["yt/utilities/lib/origami.pyx",
"yt/utilities/lib/origami_tags.c"],
@@ -283,6 +286,7 @@
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],
depends=["yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/mesh_intersection.pxd",
"yt/utlilites/lib/mesh_samplers.pxd",
"yt/utlilites/lib/mesh_traversal.pxd"]),
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -76,10 +76,15 @@
field_data = ad[field]
buff_size = size[0:dimension] + (1,) + size[dimension:]
+ ax = data_source.axis
+ xax = self.x_axis[ax]
+ yax = self.y_axis[ax]
c = np.float64(data_source.center[dimension].d)
- bounds.insert(2*dimension, c)
- bounds.insert(2*dimension, c)
- bounds = np.reshape(bounds, (3, 2))
+
+ extents = np.zeros((3, 2))
+ extents[ax] = np.array([c, c])
+ extents[xax] = bounds[0:2]
+ extents[yax] = bounds[2:4]
# if this is an element field, promote to 2D here
if len(field_data.shape) == 1:
@@ -101,13 +106,10 @@
img = pixelize_element_mesh(coords,
indices,
- buff_size, field_data, bounds,
+ buff_size, field_data, extents,
index_offset=offset)
# re-order the array and squeeze out the dummy dim
- ax = data_source.axis
- xax = self.x_axis[ax]
- yax = self.y_axis[ax]
return np.squeeze(np.transpose(img, (yax, xax, ax)))
elif dimension < 3:
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -4,7 +4,7 @@
from yt.utilities.lib.element_mappings cimport ElementSampler
from yt.utilities.lib.primitives cimport BBox, Ray
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -30,11 +30,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -52,11 +47,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef class P1Sampler3D(ElementSampler):
@@ -72,11 +62,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -169,11 +154,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -191,11 +171,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -214,11 +189,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -250,8 +220,3 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
-
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -86,15 +86,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- pass
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
pass
@@ -182,19 +173,6 @@
return 0
return 1
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
cdef class P1Sampler3D(ElementSampler):
'''
@@ -261,19 +239,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double u, v, w
cdef double thresh = 2.0e-2
@@ -413,18 +378,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
@@ -568,19 +521,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -781,22 +721,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (direction == 2 and (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance)):
- return 1
- elif (direction == 1 and (fabs(mapped_coord[direction]) < tolerance)):
- return 1
- elif (direction == 0 and (fabs(mapped_coord[0]) < tolerance or
- fabs(mapped_coord[0] + mapped_coord[1] - 1.0) < tolerance)):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double r, s, t
cdef double thresh = 5.0e-2
@@ -969,20 +893,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -18,8 +18,9 @@
cimport cython
from libc.stdlib cimport malloc, free
from yt.utilities.lib.fp_utils cimport fclip, i64clip
-from libc.math cimport copysign
+from libc.math cimport copysign, fabs
from yt.utilities.exceptions import YTDomainOverflow
+from yt.utilities.lib.vec3_ops cimport subtract, cross, dot, L2_norm
cdef extern from "math.h":
double exp(double x) nogil
@@ -501,6 +502,11 @@
cdef np.float64_t p0[3]
cdef np.float64_t p1[3]
cdef np.float64_t p3[3]
+ cdef np.float64_t E0[3]
+ cdef np.float64_t E1[3]
+ cdef np.float64_t tri_norm[3]
+ cdef np.float64_t plane_norm[3]
+ cdef np.float64_t dp
cdef int i, j, k, count, ntri, nlines
nlines = 0
ntri = triangles.shape[0]
@@ -510,6 +516,22 @@
first = last = points = NULL
for i in range(ntri):
count = 0
+
+ # skip if triangle is close to being parallel to plane
+ for j in range(3):
+ p0[j] = triangles[i, 0, j]
+ p1[j] = triangles[i, 1, j]
+ p3[j] = triangles[i, 2, j]
+ plane_norm[j] = 0.0
+ plane_norm[ax] = 1.0
+ subtract(p1, p0, E0)
+ subtract(p3, p0, E1)
+ cross(E0, E1, tri_norm)
+ dp = dot(tri_norm, plane_norm)
+ dp /= L2_norm(tri_norm)
+ if (fabs(dp) > 0.995):
+ continue
+
# Now for each line segment (01, 12, 20) we check to see how many cross
# the coordinate.
for j in range(3):
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_construction.h
--- a/yt/utilities/lib/mesh_construction.h
+++ /dev/null
@@ -1,66 +0,0 @@
-#define MAX_NUM_TRI 12
-#define HEX_NV 8
-#define HEX_NT 12
-#define TETRA_NV 4
-#define TETRA_NT 4
-#define WEDGE_NV 6
-#define WEDGE_NT 8
-
-// This array is used to triangulate the hexahedral mesh elements
-// Each element has six faces with two triangles each.
-// The vertex ordering convention is assumed to follow that used
-// here: http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-// Note that this is the case for Exodus II data.
-int triangulate_hex[MAX_NUM_TRI][3] = {
- {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
- {4, 5, 6}, {4, 6, 7}, // Face is 4 5 6 7
- {0, 1, 5}, {0, 5, 4}, // Face is 0 1 5 4
- {1, 2, 6}, {1, 6, 5}, // Face is 1 2 6 5
- {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
- {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
-};
-
-// Similarly, this is used to triangulate the tetrahedral cells
-int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 1, 3},
- {0, 2, 3},
- {1, 2, 3},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-// Triangulate wedges
-int triangulate_wedge[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 3, 1},
- {1, 3, 4},
- {0, 2, 3},
- {2, 5, 3},
- {1, 4, 2},
- {2, 4, 5},
- {3, 5, 4},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-
-// This is used to select faces from a 20-sided hex element
-int hex20_faces[6][8] = {
- {0, 1, 5, 4, 12, 8, 13, 16},
- {1, 2, 6, 5, 13, 9, 14, 17},
- {3, 2, 6, 7, 15, 10, 14, 18},
- {0, 3, 7, 4, 12, 11, 15, 19},
- {4, 5, 6, 7, 19, 16, 17, 18},
- {0, 1, 2, 3, 11, 8, 9, 10}
-};
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -40,7 +40,7 @@
patchBoundsFunc
from yt.utilities.exceptions import YTElementTypeNotRecognized
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_triangulation.h
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -0,0 +1,66 @@
+#define MAX_NUM_TRI 12
+#define HEX_NV 8
+#define HEX_NT 12
+#define TETRA_NV 4
+#define TETRA_NT 4
+#define WEDGE_NV 6
+#define WEDGE_NT 8
+
+// This array is used to triangulate the hexahedral mesh elements
+// Each element has six faces with two triangles each.
+// The vertex ordering convention is assumed to follow that used
+// here: http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
+// Note that this is the case for Exodus II data.
+int triangulate_hex[MAX_NUM_TRI][3] = {
+ {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
+ {4, 5, 6}, {4, 6, 7}, // Face is 4 5 6 7
+ {0, 1, 5}, {0, 5, 4}, // Face is 0 1 5 4
+ {1, 2, 6}, {1, 6, 5}, // Face is 1 2 6 5
+ {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
+ {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
+};
+
+// Similarly, this is used to triangulate the tetrahedral cells
+int triangulate_tetra[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 1, 3},
+ {0, 2, 3},
+ {1, 2, 3},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+// Triangulate wedges
+int triangulate_wedge[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 3, 1},
+ {1, 3, 4},
+ {0, 2, 3},
+ {2, 5, 3},
+ {1, 4, 2},
+ {2, 4, 5},
+ {3, 5, 4},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+
+// This is used to select faces from a 20-sided hex element
+int hex20_faces[6][8] = {
+ {0, 1, 5, 4, 12, 8, 13, 16},
+ {1, 2, 6, 5, 13, 9, 14, 17},
+ {3, 2, 6, 7, 15, 10, 14, 18},
+ {0, 3, 7, 4, 12, 11, 15, 19},
+ {4, 5, 6, 7, 19, 16, 17, 18},
+ {0, 1, 2, 3, 11, 8, 9, 10}
+};
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_triangulation.pyx
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -0,0 +1,58 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+
+from yt.utilities.exceptions import YTElementTypeNotRecognized
+
+cdef extern from "mesh_triangulation.h":
+ enum:
+ MAX_NUM_TRI
+
+ int HEX_NV
+ int HEX_NT
+ int TETRA_NV
+ int TETRA_NT
+ int WEDGE_NV
+ int WEDGE_NT
+ int triangulate_hex[MAX_NUM_TRI][3]
+ int triangulate_tetra[MAX_NUM_TRI][3]
+ int triangulate_wedge[MAX_NUM_TRI][3]
+ int hex20_faces[6][8]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
+ cdef np.int64_t num_elem = indices.shape[0]
+ cdef np.int64_t VPE = indices.shape[1] # num verts per element
+ cdef np.int64_t TPE # num triangles per element
+ cdef int[MAX_NUM_TRI][3] tri_array
+
+ if (VPE == 8 or VPE == 20 or VPE == 27):
+ TPE = HEX_NT
+ tri_array = triangulate_hex
+ elif VPE == 4:
+ TPE = TETRA_NT
+ tri_array = triangulate_tetra
+ elif VPE == 6:
+ TPE = WEDGE_NT
+ tri_array = triangulate_wedge
+ else:
+ raise YTElementTypeNotRecognized(3, VPE)
+
+ cdef np.int64_t num_tri = TPE * num_elem
+
+ cdef np.ndarray[np.int64_t, ndim=2] tri_indices
+ tri_indices = np.empty((num_tri, 3), dtype="int64")
+
+ cdef np.int64_t *tri_indices_ptr = <np.int64_t*> tri_indices.data
+ cdef np.int64_t *indices_ptr = <np.int64_t*> indices.data
+
+ cdef np.int64_t i, j, k, offset
+ for i in range(num_elem):
+ for j in range(TPE):
+ for k in range(3):
+ offset = tri_array[j][k]
+ tri_indices_ptr[i*TPE*3 + 3*j + k] = indices_ptr[i*VPE + offset]
+ return tri_indices
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -18,7 +18,8 @@
cimport cython
cimport libc.math as math
from yt.utilities.lib.fp_utils cimport fmin, fmax, i64min, i64max, imin, imax, fabs
-from yt.utilities.exceptions import YTPixelizeError, \
+from yt.utilities.exceptions import \
+ YTPixelizeError, \
YTElementTypeNotRecognized
from libc.stdlib cimport malloc, free
from vec3_ops cimport dot, cross, subtract
@@ -544,8 +545,7 @@
buff_size,
np.ndarray[np.float64_t, ndim=2] field,
extents,
- int index_offset = 0,
- double thresh = -1.0):
+ int index_offset = 0):
cdef np.ndarray[np.float64_t, ndim=3] img
img = np.zeros(buff_size, dtype="float64")
# Two steps:
@@ -592,27 +592,10 @@
# if we are in 2D land, the 1 cell thick dimension had better be 'z'
if ndim == 2:
- assert(buff_size[2] == 1)
+ if buff_size[2] != 1:
+ raise RuntimeError("Slices of 2D datasets must be "
+ "perpendicular to the 'z' direction.")
- ax = -1
- for i in range(3):
- if buff_size[i] == 1:
- ax = i
- if ax == -1:
- raise RuntimeError
- xax = yax = -1
- if ax == 0:
- xax = 1
- yax = 2
- elif ax == 1:
- xax = 1
- yax = 2
- elif ax == 2:
- xax = 0
- yax = 1
- if xax == -1 or yax == -1:
- raise RuntimeError
-
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
@@ -673,20 +656,11 @@
sampler.map_real_to_unit(mapped_coord, vertices, ppoint)
if not sampler.check_inside(mapped_coord):
continue
- # if thresh is negative, we do the normal sample operation
- if thresh < 0.0:
- if (num_field_vals == 1):
- img[pi, pj, pk] = field_vals[0]
- else:
- img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
- field_vals)
+ if (num_field_vals == 1):
+ img[pi, pj, pk] = field_vals[0]
else:
- # otherwise, we draw the element boundaries
- if sampler.check_near_edge(mapped_coord, thresh, xax):
- img[pi, pj, pk] = 1.0
- elif sampler.check_near_edge(mapped_coord, thresh, yax):
- img[pi, pj, pk] = 1.0
-
+ img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
+ field_vals)
free(vertices)
free(field_vals)
return img
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -22,8 +22,6 @@
from distutils.version import LooseVersion
-from yt.config import \
- ytcfg
from yt.funcs import \
mylog, iterable
from yt.extern.six import add_metaclass
@@ -31,13 +29,14 @@
from yt.visualization.image_writer import apply_colormap
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.utilities.lib.pixelization_routines import \
- pixelize_element_mesh, pixelize_off_axis_cartesian, \
+ pixelize_off_axis_cartesian, \
pixelize_cartesian
from yt.analysis_modules.cosmological_observation.light_ray.light_ray import \
periodic_ray
from yt.utilities.lib.line_integral_convolution import \
line_integral_convolution_2d
-
+from yt.geometry.unstructured_mesh_handler import UnstructuredIndex
+from yt.utilities.lib.mesh_triangulation import triangulate_indices
from . import _MPL
@@ -1605,85 +1604,77 @@
class MeshLinesCallback(PlotCallback):
"""
- annotate_mesh_lines(thresh=0.01)
+ annotate_mesh_lines()
- Adds the mesh lines to the plot.
-
- This is done by marking the pixels where the mapped coordinate
- is within some threshold distance of one of the element boundaries.
- If the mesh lines are too thick or too thin, try varying thresh.
+ Adds mesh lines to the plot. Only works for unstructured or
+ semi-structured mesh data. For structured grid data, see
+ GridBoundaryCallback or CellEdgesCallback.
Parameters
----------
- thresh : float
- The threshold distance, in mapped coordinates, within which the
- pixels will be marked as part of the element boundary.
- Default is 0.01.
+ plot_args: dict, optional
+ A dictionary of arguments that will be passed to matplotlib.
+
+ Example
+ -------
+
+ >>> import yt
+ >>> ds = yt.load("MOOSE_sample_data/out.e-s010")
+ >>> sl = yt.SlicePlot(ds, 'z', ('connect2', 'convected'))
+ >>> sl.annotate_mesh_lines(plot_args={'color':'black'})
"""
_type_name = "mesh_lines"
- def __init__(self, thresh=0.1, cmap=None):
+ def __init__(self, plot_args=None):
super(MeshLinesCallback, self).__init__()
- self.thresh = thresh
- if cmap is None:
- cmap = ytcfg.get("yt", "default_colormap")
- self.cmap = cmap
+ self.plot_args = plot_args
+
+ def promote_2d_to_3d(self, coords, indices, plot):
+ new_coords = np.zeros((2*coords.shape[0], 3))
+ new_connects = np.zeros((indices.shape[0], 2*indices.shape[1]),
+ dtype=np.int64)
+
+ new_coords[0:coords.shape[0],0:2] = coords
+ new_coords[0:coords.shape[0],2] = plot.ds.domain_left_edge[2]
+ new_coords[coords.shape[0]:,0:2] = coords
+ new_coords[coords.shape[0]:,2] = plot.ds.domain_right_edge[2]
+
+ new_connects[:,0:indices.shape[1]] = indices
+ new_connects[:,indices.shape[1]:] = indices + coords.shape[0]
+
+ return new_coords, new_connects
def __call__(self, plot):
- ftype, fname = plot.field
- mesh_id = int(ftype[-1]) - 1
- mesh = plot.ds.index.meshes[mesh_id]
+ index = plot.ds.index
+ if not issubclass(type(index), UnstructuredIndex):
+ raise RuntimeError("Mesh line annotations only work for "
+ "unstructured or semi-structured mesh data.")
+ for i, m in enumerate(index.meshes):
+ try:
+ ftype, fname = plot.field
+ if ftype.startswith('connect') and int(ftype[-1]) - 1 != i:
+ continue
+ except ValueError:
+ pass
+ coords = m.connectivity_coords
+ indices = m.connectivity_indices - m._index_offset
+
+ num_verts = indices.shape[1]
+ num_dims = coords.shape[1]
- coords = mesh.connectivity_coords
- indices = mesh.connectivity_indices
+ if num_dims == 2 and num_verts == 3:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
+ elif num_dims == 2 and num_verts == 4:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
- xx0, xx1 = plot._axes.get_xlim()
- yy0, yy1 = plot._axes.get_ylim()
-
- offset = mesh._index_offset
- ax = plot.data.axis
- xax = plot.data.ds.coordinates.x_axis[ax]
- yax = plot.data.ds.coordinates.y_axis[ax]
-
- size = plot.image._A.shape
- c = plot.data.center[ax]
- buff_size = size[0:ax] + (1,) + size[ax:]
-
- x0, x1 = plot.xlim
- y0, y1 = plot.ylim
- c = plot.data.center[ax]
- bounds = [x0, x1, y0, y1]
- bounds.insert(2*ax, c)
- bounds.insert(2*ax, c)
- bounds = np.reshape(bounds, (3, 2))
-
- # draw the mesh lines by marking where the mapped
- # coordinate is close to the domain edge. We do this by
- # calling the pixelizer with a dummy field and passing in
- # a non-negative thresh.
- dummy_field = np.zeros(indices.shape, dtype=np.float64)
- img = pixelize_element_mesh(coords, indices,
- buff_size,
- dummy_field,
- bounds,
- offset, self.thresh)
- img = np.squeeze(np.transpose(img, (yax, xax, ax)))
-
- # convert to RGBA
- size = plot.frb.buff_size
- image = np.zeros((size[0], size[1], 4), dtype=np.uint8)
- image[:, :, 0][img > 0.0] = 0
- image[:, :, 1][img > 0.0] = 0
- image[:, :, 2][img > 0.0] = 0
- image[:, :, 3][img > 0.0] = 255
-
- plot._axes.imshow(image, zorder=1,
- extent=[xx0, xx1, yy0, yy1],
- origin='lower', cmap=self.cmap,
- interpolation='nearest')
+ tri_indices = triangulate_indices(indices)
+ points = coords[tri_indices]
+
+ tfc = TriangleFacetsCallback(points, plot_args=self.plot_args)
+ tfc(plot)
class TriangleFacetsCallback(PlotCallback):
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -21,7 +21,9 @@
from yt.config import \
ytcfg
from yt.testing import \
- fake_amr_ds
+ fake_amr_ds, \
+ fake_tetrahedral_ds, \
+ fake_hexahedral_ds
import yt.units as u
from .test_plotwindow import assert_fname
from yt.utilities.exceptions import \
@@ -368,6 +370,22 @@
color=(0.0, 1.0, 1.0))
p.save(prefix)
+def test_mesh_lines_callback():
+ with _cleanup_fname() as prefix:
+
+ ds = fake_hexahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+ ds = fake_tetrahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+
def test_line_integral_convolution_callback():
with _cleanup_fname() as prefix:
ds = fake_amr_ds(fields =
@@ -390,3 +408,4 @@
alpha=0.9, const_alpha=True)
p.save(prefix)
+
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