[yt-svn] commit/yt: 20 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Oct 17 14:20:13 PDT 2016
20 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/d0f56ef727ea/
Changeset: d0f56ef727ea
Branch: yt
User: al007
Date: 2016-09-27 18:18:40+00:00
Summary: Comment out stripping to first order for second order tets.
Affected #: 2 files
diff -r 5a08f4cc1af1a18889e0dc93cb3afa908f2d3d2f -r d0f56ef727eabada948800eba814d03ffb981782 setup.py
--- a/setup.py
+++ b/setup.py
@@ -322,7 +322,7 @@
def finalize_options(self):
from Cython.Build import cythonize
self.distribution.ext_modules[:] = cythonize(
- self.distribution.ext_modules)
+ self.distribution.ext_modules, gdb_debug=True)
_build_ext.finalize_options(self)
# Prevent numpy from thinking it is still in its setup process
# see http://stackoverflow.com/a/21621493/1382869
@@ -337,7 +337,7 @@
def run(self):
# Make sure the compiled Cython files in the distribution are up-to-date
from Cython.Build import cythonize
- cythonize(cython_extensions)
+ cythonize(cython_extensions, gdb_dbg=True)
_sdist.run(self)
setup(
diff -r 5a08f4cc1af1a18889e0dc93cb3afa908f2d3d2f -r d0f56ef727eabada948800eba814d03ffb981782 yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -131,7 +131,7 @@
--------
The easiest way to make a VolumeSource is to use the volume_render
- function, so that the VolumeSource gets created automatically. This
+ function, so that the VolumeSource gets created automatically. This
example shows how to do this and then access the resulting source:
>>> import yt
@@ -545,7 +545,7 @@
'''
This is the name of the colormap that will be used when rendering
this MeshSource object. Should be a string, like 'arbre', or 'dusk'.
-
+
'''
def fget(self):
@@ -562,7 +562,7 @@
'''
These are the bounds that will be used with the colormap to the display
the rendered image. Should be a (vmin, vmax) tuple, like (0.0, 2.0). If
- None, the bounds will be automatically inferred from the max and min of
+ None, the bounds will be automatically inferred from the max and min of
the rendered data.
'''
@@ -603,10 +603,10 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry. Right now, high-order geometry is only
# implemented for 20-point hexes.
- if indices.shape[1] == 20:
+ if indices.shape[1] == 20 or indices.shape[1] == 10:
self.mesh = mesh_construction.QuadraticElementMesh(self.volume,
vertices,
indices,
@@ -620,12 +620,12 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
+ # elif indices.shape[1] == 10:
+ # # tetrahedral
+ # mylog.warning("10-node tetrahedral elements not yet supported, " +
+ # "dropping to 1st order.")
+ # field_data = field_data[:, 0:4]
+ # indices = indices[:, 0:4]
self.mesh = mesh_construction.LinearElementMesh(self.volume,
vertices,
@@ -651,7 +651,7 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry.
if indices.shape[1] == 27:
# hexahedral
@@ -659,12 +659,12 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
+ # elif indices.shape[1] == 10:
+ # # tetrahedral
+ # mylog.warning("10-node tetrahedral elements not yet supported, " +
+ # "dropping to 1st order.")
+ # field_data = field_data[:, 0:4]
+ # indices = indices[:, 0:4]
self.volume = BVH(vertices, indices, field_data)
@@ -794,7 +794,7 @@
cmap_name=self._cmap)/255.
alpha = image[:, :, 3]
alpha[self.sampler.aimage_used == -1] = 0.0
- image[:, :, 3] = alpha
+ image[:, :, 3] = alpha
return image
def __repr__(self):
@@ -854,7 +854,7 @@
assert(positions.ndim == 2 and positions.shape[1] == 3)
if colors is not None:
assert(colors.ndim == 2 and colors.shape[1] == 4)
- assert(colors.shape[0] == positions.shape[0])
+ assert(colors.shape[0] == positions.shape[0])
self.positions = positions
# If colors aren't individually set, make black with full opacity
if colors is None:
@@ -1057,7 +1057,7 @@
Examples
--------
- This example shows how to use BoxSource to add an outline of the
+ This example shows how to use BoxSource to add an outline of the
domain boundaries to a volume rendering.
>>> import yt
@@ -1065,12 +1065,12 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> box_source = BoxSource(ds.domain_left_edge,
... ds.domain_right_edge,
... [1.0, 1.0, 1.0, 1.0])
>>> sc.add_source(box_source)
- >>>
+ >>>
>>> im = sc.render()
"""
@@ -1078,7 +1078,7 @@
assert(left_edge.shape == (3,))
assert(right_edge.shape == (3,))
-
+
if color is None:
color = np.array([1.0, 1.0, 1.0, 1.0])
@@ -1118,7 +1118,7 @@
Examples
--------
- This example makes a volume rendering and adds outlines of all the
+ This example makes a volume rendering and adds outlines of all the
AMR grids in the simulation:
>>> import yt
@@ -1140,12 +1140,12 @@
>>> import yt
>>> from yt.visualization.volume_rendering.api import GridSource
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> dd = ds.sphere("c", (0.1, "unitary"))
>>> grid_source = GridSource(dd, alpha=1.0)
- >>>
+ >>>
>>> sc.add_source(grid_source)
>>>
>>> im = sc.render()
@@ -1223,11 +1223,11 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> coord_source = CoordinateVectorSource()
- >>>
+ >>>
>>> sc.add_source(coord_source)
- >>>
+ >>>
>>> im = sc.render()
"""
https://bitbucket.org/yt_analysis/yt/commits/d9051c62a497/
Changeset: d9051c62a497
Branch: yt
User: al007
Date: 2016-09-27 20:24:45+00:00
Summary: Prevent segmentation fault when NotImplemented error raised in bounding_volume_hierarchy
Affected #: 1 file
diff -r d0f56ef727eabada948800eba814d03ffb981782 -r d9051c62a497880c5c5915c9ab86006d0d72e458 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -45,11 +45,11 @@
'''
This class implements a bounding volume hierarchy (BVH), a spatial acceleration
- structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
- partitioning the *volume* of the parent to create the children, we partition the
+ structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
+ partitioning the *volume* of the parent to create the children, we partition the
triangles themselves into 'left' or 'right' sub-trees. The bounding volume for a
node is then determined by computing the bounding volume of the triangles that
- belong to it. This allows us to quickly discard triangles that are not close
+ belong to it. This allows us to quickly discard triangles that are not close
to intersecting a given ray.
'''
@@ -109,7 +109,7 @@
self.vertices[vertex_offset + k] = vertices[indices[i,j]][k]
field_offset = i*self.num_field_per_elem
for j in range(self.num_field_per_elem):
- self.field_data[field_offset + j] = field_data[i][j]
+ self.field_data[field_offset + j] = field_data[i][j]
# set up primitives
if self.num_verts_per_elem == 20:
@@ -124,7 +124,7 @@
self.get_bbox = triangle_bbox
self.get_intersect = ray_triangle_intersect
self._set_up_triangles(vertices, indices)
-
+
self.root = self._recursive_build(0, self.num_prim)
@cython.boundscheck(False)
@@ -181,7 +181,7 @@
tri_index,
self.centroids[tri_index])
self.get_bbox(self.primitives,
- tri_index,
+ tri_index,
&(self.bboxes[tri_index]))
cdef void _recursive_free(self, BVHNode* node) nogil:
@@ -191,6 +191,8 @@
free(node)
def __dealloc__(self):
+ if self.root == NULL:
+ return
self._recursive_free(self.root)
free(self.primitives)
free(self.prim_ids)
@@ -224,11 +226,11 @@
mid += 1
begin += 1
return mid
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil:
cdef np.int64_t i, j
cdef BBox box = self.bboxes[begin]
@@ -236,7 +238,7 @@
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],
self.bboxes[i].left_edge[j])
- box.right_edge[j] = fmax(box.right_edge[j],
+ box.right_edge[j] = fmax(box.right_edge[j],
self.bboxes[i].right_edge[j])
node.bbox = box
@@ -245,7 +247,7 @@
@cython.cdivision(True)
cdef void intersect(self, Ray* ray) nogil:
self._recursive_intersect(ray, self.root)
-
+
if ray.elem_id < 0:
return
@@ -253,7 +255,7 @@
cdef np.int64_t i
for i in range(3):
position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
-
+
cdef np.float64_t* vertex_ptr
cdef np.float64_t* field_ptr
vertex_ptr = self.vertices + ray.elem_id*self.num_verts_per_elem*3
@@ -297,17 +299,17 @@
node.end = end
self._get_node_bbox(node, begin, end)
-
+
# check for leaf
if (end - begin) <= LEAF_SIZE:
return node
-
+
# we use the "split in the middle of the longest axis approach"
# see: http://www.vadimkravcenko.com/bvh-tree-building/
# compute longest dimension
cdef np.int64_t ax = 0
- cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
+ cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
node.bbox.left_edge[0])
if fabs(node.bbox.right_edge[1] - node.bbox.left_edge[1]) > d:
ax = 1
@@ -323,7 +325,7 @@
if(mid == begin or mid == end):
mid = begin + (end-begin)/2
-
+
# recursively build sub-trees
node.left = self._recursive_build(begin, mid)
node.right = self._recursive_build(mid, end)
@@ -337,16 +339,16 @@
cdef void cast_rays(np.float64_t* image,
const np.float64_t* origins,
const np.float64_t* direction,
- const int N,
+ const int N,
BVH bvh) nogil:
- cdef Ray* ray
+ cdef Ray* ray
cdef int i, j, k
-
+
with nogil, parallel():
-
+
ray = <Ray *> malloc(sizeof(Ray))
-
+
for k in range(3):
ray.direction[k] = direction[k]
ray.inv_dir[k] = 1.0 / direction[k]
@@ -366,11 +368,11 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
+def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
np.ndarray[np.float64_t, ndim=2] origins,
np.ndarray[np.float64_t, ndim=1] direction,
BVH bvh):
-
+
cdef int N = origins.shape[0]
cast_rays(&image[0], &origins[0, 0], &direction[0], N, bvh)
https://bitbucket.org/yt_analysis/yt/commits/d2e1c551b36c/
Changeset: d2e1c551b36c
Branch: yt
User: al007
Date: 2016-09-27 20:51:12+00:00
Summary: First crack at rendering second order tets. Triangulation probably wrong.
Affected #: 1 file
diff -r d9051c62a497880c5c5915c9ab86006d0d72e458 -r d2e1c551b36ca16daf80b3cce2a2c10c6b4cc8ce yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -24,13 +24,16 @@
Q1Sampler3D, \
P1Sampler3D, \
W1Sampler3D, \
- S2Sampler3D
+ S2Sampler3D, \
+ Tet2Sampler3D
+
from yt.utilities.lib.vec3_ops cimport L2_norm
cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
@@ -82,6 +85,10 @@
elif self.num_verts_per_elem == 20:
self.num_prim_per_elem = 6
self.sampler = S2Sampler
+ elif self.num_verts_per_elem == 10:
+ self.num_prim_per_elem = 4
+ self.tri_array = triangulate_tetra
+ self.sampler = Tet2Sampler
else:
raise NotImplementedError("Could not determine element type for "
"nverts = %d. " % self.num_verts_per_elem)
https://bitbucket.org/yt_analysis/yt/commits/c7497ba2296a/
Changeset: c7497ba2296a
Branch: yt
User: al007
Date: 2016-09-28 23:05:36+00:00
Summary: Set num_prim_per_elem to TETRA_NT
Affected #: 1 file
diff -r d2e1c551b36ca16daf80b3cce2a2c10c6b4cc8ce -r c7497ba2296aa327fc6b5c693a15307b5471f811 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -86,7 +86,7 @@
self.num_prim_per_elem = 6
self.sampler = S2Sampler
elif self.num_verts_per_elem == 10:
- self.num_prim_per_elem = 4
+ self.num_prim_per_elem = TETRA_NT
self.tri_array = triangulate_tetra
self.sampler = Tet2Sampler
else:
https://bitbucket.org/yt_analysis/yt/commits/df0fe5963732/
Changeset: df0fe5963732
Branch: yt
User: al007
Date: 2016-09-30 20:54:48+00:00
Summary: Delete redundant and troublesome switch_orientation call in position's fget method.
Affected #: 1 file
diff -r 1d70bf4124996d584483c173e37f29670a73db99 -r df0fe5963732f3e2b3e899c7a25e9cdebeb30296 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -158,7 +158,7 @@
def position():
doc = '''
- The location of the camera.
+ The location of the camera.
Parameters
----------
@@ -178,8 +178,6 @@
raise RuntimeError(
'Cannot set the camera focus and position to the same value')
self._position = position
- self.switch_orientation(normal_vector=self.focus - self._position,
- north_vector=None)
def fdel(self):
del self._position
@@ -483,11 +481,11 @@
>>> sc = Scene()
>>> cam = sc.add_camera()
>>> # rotate the camera by pi / 4 radians:
- >>> cam.rotate(np.pi/4.0)
+ >>> cam.rotate(np.pi/4.0)
>>> # rotate the camera about the y-axis instead of cam.north_vector:
- >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
>>> # rotate the camera about the origin instead of its own position:
- >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
"""
rotate_all = rot_vector is not None
@@ -593,7 +591,7 @@
>>> sc = Scene()
>>> cam = sc.add_camera(ds)
>>> # roll the camera by pi / 4 radians:
- >>> cam.roll(np.pi/4.0)
+ >>> cam.roll(np.pi/4.0)
>>> # roll the camera about the origin instead of its own position:
>>> cam.roll(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
@@ -627,7 +625,7 @@
>>> import yt
>>> import numpy as np
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
>>> cam = sc.camera
>>> for i in cam.iter_rotate(np.pi, 10):
@@ -688,7 +686,7 @@
def zoom(self, factor):
r"""Change the width of the FOV of the camera.
- This will appear to zoom the camera in by some `factor` toward the
+ This will appear to zoom the camera in by some `factor` toward the
focal point along the current view direction, but really it's just
changing the width of the field of view.
https://bitbucket.org/yt_analysis/yt/commits/c396959f0d97/
Changeset: c396959f0d97
Branch: yt
User: al007
Date: 2016-09-30 21:06:28+00:00
Summary: Add cam_pos and north_vector test. Closes #1287
Affected #: 1 file
diff -r df0fe5963732f3e2b3e899c7a25e9cdebeb30296 -r c396959f0d975a327974d8e2667aa000acd2c8d5 yt/visualization/volume_rendering/tests/test_camera_attributes.py
--- a/yt/visualization/volume_rendering/tests/test_camera_attributes.py
+++ b/yt/visualization/volume_rendering/tests/test_camera_attributes.py
@@ -110,3 +110,11 @@
for lens_type in valid_lens_types:
cam.set_lens(lens_type)
+
+ cam.focus = [0, 0, 0]
+ cam_pos = [1, 0, 0]
+ north_vector = [0, 1, 0]
+ cam.set_position(cam_pos, north_vector)
+ cam_pos = [0, 1, 0]
+ north_vector = [0, 0, 1]
+ cam.set_position(cam_pos, north_vector)
https://bitbucket.org/yt_analysis/yt/commits/2164544a5466/
Changeset: 2164544a5466
Branch: yt
User: al007
Date: 2016-09-30 21:24:59+00:00
Summary: Swap which switch_orientation call is removed.
Affected #: 1 file
diff -r c396959f0d975a327974d8e2667aa000acd2c8d5 -r 2164544a546653def49e60a18c0d5238d90ff716 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -178,6 +178,8 @@
raise RuntimeError(
'Cannot set the camera focus and position to the same value')
self._position = position
+ self.switch_orientation(normal_vector=self.focus - self._position,
+ north_vector=self.north_vector)
def fdel(self):
del self._position
@@ -384,9 +386,9 @@
calculated automatically.
"""
+ if north_vector is not None:
+ self.north_vector = north_vector
self.position = position
- self.switch_orientation(normal_vector=self.focus - self.position,
- north_vector=north_vector)
def get_position(self):
"""Return the current camera position"""
https://bitbucket.org/yt_analysis/yt/commits/38390ecb88cf/
Changeset: 38390ecb88cf
Branch: yt
User: al007
Date: 2016-09-30 21:27:21+00:00
Summary: Added issue reference comment to test.
Affected #: 1 file
diff -r 2164544a546653def49e60a18c0d5238d90ff716 -r 38390ecb88cf22ba7ffadfb033626913b14e9273 yt/visualization/volume_rendering/tests/test_camera_attributes.py
--- a/yt/visualization/volume_rendering/tests/test_camera_attributes.py
+++ b/yt/visualization/volume_rendering/tests/test_camera_attributes.py
@@ -111,6 +111,7 @@
for lens_type in valid_lens_types:
cam.set_lens(lens_type)
+ # See issue #1287
cam.focus = [0, 0, 0]
cam_pos = [1, 0, 0]
north_vector = [0, 1, 0]
https://bitbucket.org/yt_analysis/yt/commits/05262ce49422/
Changeset: 05262ce49422
Branch: yt
User: al007
Date: 2016-09-30 16:35:46+00:00
Summary: Adding tet_patches.
Affected #: 7 files
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
# node for the bounding volume hierarchy
cdef struct BVHNode:
@@ -70,7 +71,7 @@
np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil
cdef void intersect(self, Ray* ray) nogil
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil
cdef void _recursive_intersect(self, Ray* ray, BVHNode* node) nogil
cdef BVHNode* _recursive_build(self, np.int64_t begin, np.int64_t end) nogil
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -18,7 +18,12 @@
Patch, \
ray_patch_intersect, \
patch_centroid, \
- patch_bbox
+ patch_bbox, \
+ TetPatch, \
+ ray_tetPatch_intersect, \
+ tetPatch_centroid, \
+ tetPatch_bbox
+
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
@@ -86,8 +91,7 @@
self.num_prim_per_elem = 6
self.sampler = S2Sampler
elif self.num_verts_per_elem == 10:
- self.num_prim_per_elem = TETRA_NT
- self.tri_array = triangulate_tetra
+ self.num_prim_per_elem = 4
self.sampler = Tet2Sampler
else:
raise NotImplementedError("Could not determine element type for "
@@ -125,6 +129,12 @@
self.get_bbox = patch_bbox
self.get_intersect = ray_patch_intersect
self._set_up_patches(vertices, indices)
+ elif self.num_verts_per_elem == 10:
+ self.primitives = malloc(self.num_prim * sizeof(TetPatch))
+ self.get_centroid = tetPatch_centroid
+ self.get_bbox = tetPatch_bbox
+ self.get_intersect = ray_tetPatch_intersect
+ self._set_up_tetPatches(vertices, indices)
else:
self.primitives = malloc(self.num_prim * sizeof(Triangle))
self.get_centroid = triangle_centroid
@@ -163,6 +173,32 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+ cdef void _set_up_tetPatches(self, np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil:
+ cdef TetPatch* tetPatch
+ cdef np.int64_t i, j, k, ind, idim
+ cdef np.int64_t offset, prim_index
+ for i in range(self.num_elem):
+ offset = self.num_prim_per_elem*i
+ for j in range(self.num_prim_per_elem): # for each face
+ prim_index = offset + j
+ tetPatch = &( <TetPatch*> self.primitives)[prim_index]
+ self.prim_ids[prim_index] = prim_index
+ tetPatch.elem_id = i
+ for k in range(6): # for each vertex
+ ind = tet10_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ tetPatch.v[k][idim] = vertices[indices[i, ind]][idim]
+ self.get_centroid(self.primitives,
+ prim_index,
+ self.centroids[prim_index])
+ self.get_bbox(self.primitives,
+ prim_index,
+ &(self.bboxes[prim_index]))
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
cdef void _set_up_triangles(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:
# fill our array of primitives
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -1,5 +1,5 @@
"""
-This file contains the ElementMesh classes, which represent the target that the
+This file contains the ElementMesh classes, which represent the target that the
rays will be cast at when rendering finite element data. This class handles
the interface between the internal representation of the mesh and the pyembree
representation.
@@ -21,7 +21,7 @@
from libc.math cimport fmax, sqrt
cimport numpy as np
-cimport pyembree.rtcore as rtc
+cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu
@@ -43,7 +43,7 @@
cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
-
+
int HEX_NV
int HEX_NT
int TETRA_NV
@@ -60,7 +60,7 @@
r'''
This creates a 1st-order mesh to be ray-traced with embree.
- Currently, we handle non-triangular mesh types by converting them
+ Currently, we handle non-triangular mesh types by converting them
to triangular meshes. This class performs this transformation.
Currently, this is implemented for hexahedral and tetrahedral
meshes.
@@ -71,21 +71,21 @@
scene : EmbreeScene
This is the scene to which the constructed polygons will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the polygon mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the polygon mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
- This should either have the shape (num_elements, 4) or
- (num_elements, 8) for tetrahedral and hexahedral meshes,
- respectively. For tetrahedral meshes, each element will
+ This should either have the shape (num_elements, 4) or
+ (num_elements, 8) for tetrahedral and hexahedral meshes,
+ respectively. For tetrahedral meshes, each element will
be represented by four triangles in the scene. For hex meshes,
- each element will be represented by 12 triangles, 2 for each
+ each element will be represented by 12 triangles, 2 for each
face. For hex meshes, we assume that the node ordering is as
- defined here:
+ defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
cdef Vertex* vertices
@@ -93,7 +93,7 @@
cdef unsigned int mesh
cdef double* field_data
cdef rtcg.RTCFilterFunc filter_func
- # triangles per element, vertices per element, and field points per
+ # triangles per element, vertices per element, and field points per
# element, respectively:
cdef int tpe, vpe, fpe
cdef int[MAX_NUM_TRI][3] tri_array
@@ -101,7 +101,7 @@
cdef MeshDataContainer datac
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray data):
@@ -119,7 +119,7 @@
self.tpe = TETRA_NT
self.tri_array = triangulate_tetra
else:
- raise YTElementTypeNotRecognized(vertices.shape[1],
+ raise YTElementTypeNotRecognized(vertices.shape[1],
indices.shape[1])
self._build_from_indices(scene, vertices, indices)
@@ -142,7 +142,7 @@
for i in range(nv):
vertices[i].x = vertices_in[i, 0]
vertices[i].y = vertices_in[i, 1]
- vertices[i].z = vertices_in[i, 2]
+ vertices[i].z = vertices_in[i, 2]
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_VERTEX_BUFFER,
vertices, 0, sizeof(Vertex))
@@ -156,7 +156,7 @@
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_INDEX_BUFFER,
triangles, 0, sizeof(Triangle))
- cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
+ cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
for i in range(ne):
for j in range(self.vpe):
element_indices[i*self.vpe + j] = indices_in[i][j]
@@ -189,7 +189,7 @@
datac.vpe = self.vpe
datac.fpe = self.fpe
self.datac = datac
-
+
rtcg.rtcSetUserData(scene.scene_i, self.mesh, &self.datac)
cdef void _set_sampler_type(self, YTEmbreeScene scene):
@@ -206,7 +206,7 @@
rtcg.rtcSetIntersectionFilterFunction(scene.scene_i,
self.mesh,
self.filter_func)
-
+
def __dealloc__(self):
free(self.field_data)
free(self.element_indices)
@@ -227,35 +227,42 @@
scene : EmbreeScene
This is the scene to which the constructed patches will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
This should have the shape (num_elements, 20). Each hex will be
- represented in the scene by 6 bi-quadratic patches. We assume that
- the node ordering is as defined here:
+ represented in the scene by 6 bi-quadratic patches. We assume that
+ the node ordering is as defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
cdef Patch* patches
cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef unsigned int mesh
- # patches per element, vertices per element, and field points per
- # element, respectively:
- cdef int ppe, vpe, fpe
+ # patches per element, vertices per element, vertices per face,
+ # and field points per element, respectively:
+ cdef int ppe, vpe, vpf, fpe
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray field_data):
- # only 20-point hexes are supported right now.
+ # 20-point hexes
if indices.shape[1] == 20:
self.vpe = 20
+ self.ppe = 6
+ self.vpf = 8
+ # 10-point tets
+ elif indices.shape[1] == 10:
+ self.vpe = 10
+ self.ppe = 4
+ self.vpf = 6
else:
raise NotImplementedError
@@ -268,33 +275,41 @@
cdef int i, j, ind, idim
cdef int ne = indices_in.shape[0]
cdef int nv = vertices_in.shape[0]
- cdef int npatch = 6*ne;
+ cdef int npatch = self.ppe*ne;
cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
cdef np.ndarray[np.float64_t, ndim=2] element_vertices
cdef Patch* patches = <Patch*> malloc(npatch * sizeof(Patch))
- self.vertices = <np.float64_t*> malloc(20 * ne * 3 * sizeof(np.float64_t))
- self.field_data = <np.float64_t*> malloc(20 * ne * sizeof(np.float64_t))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
+
+ cdef int faces[self.ppe][self.vpf];
+ if self.vpe == 20:
+ faces = hex20_faces
+ elif self.vpe == 10:
+ faces = tet10_faces
+ else:
+ raise NotImplementedError
for i in range(ne):
element_vertices = vertices_in[indices_in[i]]
- for j in range(20):
- self.field_data[i*20 + j] = field_data[i][j]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
for k in range(3):
- self.vertices[i*20*3 + j*3 + k] = element_vertices[j][k]
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
cdef Patch* patch
for i in range(ne): # for each element
element_vertices = vertices_in[indices_in[i]]
- for j in range(6): # for each face
- patch = &(patches[i*6+j])
+ for j in range(self.ppe): # for each face
+ patch = &(patches[i*self.ppe+j])
patch.geomID = mesh
- for k in range(8): # for each vertex
- ind = hex20_faces[j][k]
+ for k in range(self.vpf): # for each vertex
+ ind = faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
patch.v[k][idim] = element_vertices[ind][idim]
- patch.vertices = self.vertices + i*20*3
- patch.field_data = self.field_data + i*20
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
self.patches = patches
self.mesh = mesh
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/mesh_triangulation.h
--- a/yt/utilities/lib/mesh_triangulation.h
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -12,7 +12,7 @@
// 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
+ {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
@@ -22,7 +22,7 @@
// Similarly, this is used to triangulate the tetrahedral cells
int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 3},
+ {0, 1, 3},
{2, 3, 1},
{0, 3, 2},
{0, 2, 1},
@@ -57,10 +57,18 @@
// 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},
+ {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}
};
+
+// This is used to select faces from a second-order tet element
+int tet10_faces[4][6] = {
+ {0, 1, 3, 7, 4, 8},
+ {2, 3, 1, 5, 9, 8},
+ {0, 3, 2, 6, 7, 9},
+ {0, 2, 1, 4, 6, 5}
+};
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/mesh_triangulation.pyx
--- a/yt/utilities/lib/mesh_triangulation.pyx
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef struct TriNode:
np.uint64_t key
@@ -35,7 +36,7 @@
if not found:
return 0
return 1
-
+
cdef np.uint64_t triangle_hash(np.int64_t tri[3]) nogil:
# http://stackoverflow.com/questions/1536393/good-hash-function-for-permutations
cdef np.uint64_t h = 1
@@ -58,13 +59,13 @@
cdef TriNode **table
cdef np.uint64_t num_items
-
+
def __cinit__(self):
self.table = <TriNode**> malloc(TABLE_SIZE * sizeof(TriNode*))
for i in range(TABLE_SIZE):
self.table[i] = NULL
self.num_items = 0
-
+
def __dealloc__(self):
cdef np.int64_t i
cdef TriNode *node
@@ -77,7 +78,7 @@
free(delete_node)
self.table[i] = NULL
free(self.table)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_exterior_tris(self):
@@ -102,7 +103,7 @@
element_map[counter] = node.elem
counter += 1
node = node.next_node
-
+
return tri_indices, element_map
cdef TriNode* _allocate_new_node(self,
@@ -118,13 +119,13 @@
new_node.next_node = NULL
self.num_items += 1
return new_node
-
+
@cython.cdivision(True)
cdef void update(self, np.int64_t tri[3], np.int64_t elem) nogil:
cdef np.uint64_t key = triangle_hash(tri)
cdef np.uint64_t index = key % TABLE_SIZE
cdef TriNode *node = self.table[index]
-
+
if node == NULL:
self.table[index] = self._allocate_new_node(tri, key, elem)
return
@@ -139,7 +140,7 @@
elif node.next_node == NULL:
node.next_node = self._allocate_new_node(tri, key, elem)
return
-
+
# walk through node list
cdef TriNode* prev = node
node = node.next_node
@@ -165,7 +166,7 @@
cdef np.int64_t VPE # num verts per element
cdef np.int64_t TPE # num tris per element
cdef int[MAX_NUM_TRI][3] tri_array
-
+
def __cinit__(self, np.int64_t[:, ::1] indices):
'''
@@ -190,7 +191,7 @@
self.num_tri = self.TPE * self.num_elem
self.num_verts = self.num_tri * 3
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def cull_interior_triangles(np.int64_t[:, ::1] indices):
@@ -213,7 +214,7 @@
s.update(tri, i)
return s.get_exterior_tris()
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_vertex_data(np.float64_t[:, ::1] coords,
@@ -230,7 +231,7 @@
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t num_verts = coords.shape[0]
cdef np.float32_t[:] vertex_data = np.zeros(num_verts, dtype="float32")
-
+
cdef np.int64_t i, j
for i in range(m.num_elem):
for j in range(m.VPE):
@@ -256,17 +257,17 @@
cdef np.int64_t num_tri = exterior_tris.shape[0]
cdef np.int64_t num_verts = 3 * num_tri
cdef np.int64_t num_coords = 3 * num_verts
-
+
cdef np.float32_t[:] vertex_data
if data.ndim == 2:
vertex_data = get_vertex_data(coords, data, indices)
else:
vertex_data = data.astype("float32")
-
+
cdef np.int32_t[:] tri_indices = np.empty(num_verts, dtype=np.int32)
cdef np.float32_t[:] tri_data = np.empty(num_verts, dtype=np.float32)
cdef np.float32_t[:] tri_coords = np.empty(num_coords, dtype=np.float32)
-
+
cdef np.int64_t vert_index, i, j, k
for i in range(num_tri):
for j in range(3):
@@ -278,7 +279,7 @@
tri_indices[vert_index] = vert_index
for k in range(3):
tri_coords[vert_index*3 + k] = coords[exterior_tris[i, j], k]
-
+
return np.array(tri_coords), np.array(tri_data), np.array(tri_indices)
@cython.boundscheck(False)
@@ -291,10 +292,10 @@
coordinates and the data values.
'''
-
+
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t[:, ::1] tri_indices = np.empty((m.num_tri, 3), dtype=np.int64)
-
+
cdef np.int64_t i, j, k
for i in range(m.num_elem):
for j in range(m.TPE):
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/primitives.pxd
--- a/yt/utilities/lib/primitives.pxd
+++ b/yt/utilities/lib/primitives.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
cimport cython.floating
import numpy as np
cimport numpy as np
@@ -66,7 +66,7 @@
cdef RayHitData compute_patch_hit(cython.floating[8][3] verts,
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil
-
+
cdef np.int64_t ray_patch_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil
@@ -78,3 +78,38 @@
cdef void patch_bbox(const void *primitives,
const np.int64_t item,
BBox* bbox) nogil
+
+cdef struct TetPatch:
+ np.float64_t[6][3] v # 6 vertices per patch
+ np.int64_t elem_id
+
+cdef RayHitData compute_tetPatch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil
+
+cdef void tetPatchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil
+
+cdef void tetPatchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil
+
+cdef void tetPatchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil
+
+cdef np.int64_t ray_tetPatch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+cdef void tetPatch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+cdef void tetPatch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
diff -r c7497ba2296aa327fc6b5c693a15307b5471f811 -r 05262ce494229229ad11ec309695f445c08f0b8c yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -21,15 +21,15 @@
cdef np.float64_t tmin = -INF
cdef np.float64_t tmax = INF
-
+
cdef np.int64_t i
cdef np.float64_t t1, t2
for i in range(3):
t1 = (bbox.left_edge[i] - ray.origin[i])*ray.inv_dir[i]
- t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
+ t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
tmin = fmax(tmin, fmin(t1, t2))
tmax = fmin(tmax, fmax(t1, t2))
-
+
return tmax >= fmax(tmin, 0.0)
@cython.boundscheck(False)
@@ -53,7 +53,7 @@
cdef np.float64_t det, inv_det
det = dot(e1, P)
- if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
+ if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
return False
inv_det = 1.0 / det
@@ -87,7 +87,7 @@
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
for i in range(3):
@@ -112,7 +112,7 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceFunc(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] S) nogil:
@@ -132,9 +132,9 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivU(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
- cython.floating[3] Su) nogil:
+ cython.floating[3] Su) nogil:
cdef int i
for i in range(3):
Su[i] = (-0.25*(v - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
@@ -149,10 +149,10 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivV(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil:
- cdef int i
+ cdef int i
for i in range(3):
Sv[i] = (-0.25*(u - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
(-0.25*(u + 1.0)*(u - v - 1) + 0.25*(u + 1.0)*(v - 1.0))*verts[1][i] + \
@@ -196,7 +196,7 @@
cdef cython.floating fu = dot(N1, S) + d1
cdef cython.floating fv = dot(N2, S) + d2
cdef cython.floating err = fmax(fabs(fu), fabs(fv))
-
+
# begin Newton interation
cdef cython.floating tol = 1.0e-5
cdef int iterations = 0
@@ -213,11 +213,11 @@
J21 = dot(N2, Su)
J22 = dot(N2, Sv)
det = (J11*J22 - J12*J21)
-
+
# update the u, v values
u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
-
+
patchSurfaceFunc(verts, u, v, S)
fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
@@ -227,7 +227,7 @@
# t is the distance along the ray to this hit
cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
-
+
# return hit data
cdef RayHitData hd
hd.u = u
@@ -247,7 +247,7 @@
cdef Patch patch = (<Patch*> primitives)[item]
cdef RayHitData hd = compute_patch_hit(patch.v, ray.origin, ray.direction)
-
+
# only count this is it's the closest hit
if (hd.t < ray.t_near or hd.t > ray.t_far):
return False
@@ -259,7 +259,7 @@
return True
return False
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -291,7 +291,7 @@
cdef np.int64_t i, j
cdef Patch patch = (<Patch*> primitives)[item]
-
+
for j in range(3):
bbox.left_edge[j] = patch.v[0][j]
bbox.right_edge[j] = patch.v[0][j]
@@ -300,3 +300,196 @@
for j in range(3):
bbox.left_edge[j] = fmin(bbox.left_edge[j], patch.v[i][j])
bbox.right_edge[j] = fmax(bbox.right_edge[j], patch.v[i][j])
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tetPatchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ S[i] = (1.0 - 3.0*u + 2.0*u*u - 3.0*v + 2.0*v*v + 4.0*u*v)*verts[0][i] + \
+ (-u + 2.0*u*u)*verts[1][i] + \
+ (-v + 2.0*v*v)*verts[2][i] + \
+ (4.0*u - 4.0*u*u - 4.0*u*v)*verts[3][i] + \
+ (4.0*u*v)*verts[4][i] + \
+ (4.0*v - 4.0*v*v - 4.0*u*v)*verts[5][i]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tetPatchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil:
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Su[i] = (-3.0 + 4.0*u + 4.0*v)*verts[0][i] + \
+ (-1.0 + 4.0*u)*verts[1][i] + \
+ (4.0 - 8.0*u - 4.0*v)*verts[3][i] + \
+ (4.0*v)*verts[4][i] + \
+ (-4.0*v)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tetPatchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Sv[i] = (-3.0 + 4.0*v + 4.0*u)*verts[0][i] + \
+ (-1.0 + 4.0*v)*verts[2][i] + \
+ (-4.0*u)*verts[3][i] + \
+ (4.0*u)*verts[4][i] + \
+ (4.0 - 8.0*v - 4.0*u)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef RayHitData compute_tetPatch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil:
+
+ # first we compute the two planes that define the ray.
+ cdef cython.floating[3] n, N1, N2
+ cdef cython.floating A = dot(ray_direction, ray_direction)
+ for i in range(3):
+ n[i] = ray_direction[i] / A
+
+ if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):
+ N1[0] = n[1]
+ N1[1] =-n[0]
+ N1[2] = 0.0
+ else:
+ N1[0] = 0.0
+ N1[1] = n[2]
+ N1[2] =-n[1]
+ cross(N1, n, N2)
+
+ cdef cython.floating d1 = -dot(N1, ray_origin)
+ cdef cython.floating d2 = -dot(N2, ray_origin)
+
+ # the initial guess is set to zero
+ cdef cython.floating u = 0.0
+ cdef cython.floating v = 0.0
+ cdef cython.floating[3] S
+ tetPatchSurfaceFunc(verts, u, v, S)
+ cdef cython.floating fu = dot(N1, S) + d1
+ cdef cython.floating fv = dot(N2, S) + d2
+ cdef cython.floating err = fmax(fabs(fu), fabs(fv))
+
+ # begin Newton interation
+ cdef cython.floating tol = 1.0e-5
+ cdef int iterations = 0
+ cdef int max_iter = 10
+ cdef cython.floating[3] Su
+ cdef cython.floating[3] Sv
+ cdef cython.floating J11, J12, J21, J22, det
+ while ((err > tol) and (iterations < max_iter)):
+ # compute the Jacobian
+ tetPatchSurfaceDerivU(verts, u, v, Su)
+ tetPatchSurfaceDerivV(verts, u, v, Sv)
+ J11 = dot(N1, Su)
+ J12 = dot(N1, Sv)
+ J21 = dot(N2, Su)
+ J22 = dot(N2, Sv)
+ det = (J11*J22 - J12*J21)
+
+ # update the u, v values
+ u -= ( J22*fu - J12*fv) / det
+ v -= (-J21*fu + J11*fv) / det
+
+ tetPatchSurfaceFunc(verts, u, v, S)
+ fu = dot(N1, S) + d1
+ fv = dot(N2, S) + d2
+
+ err = fmax(fabs(fu), fabs(fv))
+ iterations += 1
+
+ # t is the distance along the ray to this hit
+ cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
+
+ # return hit data
+ cdef RayHitData hd
+ hd.u = u
+ hd.v = v
+ hd.t = t
+ hd.converged = (iterations < max_iter)
+ return hd
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_tetPatch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
+
+ cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+
+ cdef RayHitData hd = compute_tetPatch_hit(tetPatch.v, ray.origin, ray.direction)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.t_near or hd.t > ray.t_far):
+ return False
+
+ if ((0 <= hd.u <= 1.0) and (0 <= hd.v <= 1.0) and hd.converged):
+ # we have a hit, so update ray information
+ ray.t_far = hd.t
+ ray.elem_id = tetPatch.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tetPatch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ centroid[j] = 0.0
+
+ for i in range(6):
+ for j in range(3):
+ centroid[j] += tetPatch.v[i][j]
+
+ for j in range(3):
+ centroid[j] /= 6.0
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tetPatch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ bbox.left_edge[j] = tetPatch.v[0][j]
+ bbox.right_edge[j] = tetPatch.v[0][j]
+
+ for i in range(1, 6):
+ for j in range(3):
+ bbox.left_edge[j] = fmin(bbox.left_edge[j], tetPatch.v[i][j])
+ bbox.right_edge[j] = fmax(bbox.right_edge[j], tetPatch.v[i][j])
https://bitbucket.org/yt_analysis/yt/commits/4f3f870ced4e/
Changeset: 4f3f870ced4e
Branch: yt
User: al007
Date: 2016-09-30 20:49:31+00:00
Summary: Replace tetPatch with tet_patch.
Affected #: 7 files
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -70,6 +70,9 @@
cdef void _set_up_patches(self,
np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil
+ cdef void _set_up_tet_patches(self,
+ np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil
cdef void intersect(self, Ray* ray) nogil
cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -20,9 +20,9 @@
patch_centroid, \
patch_bbox, \
TetPatch, \
- ray_tetPatch_intersect, \
- tetPatch_centroid, \
- tetPatch_bbox
+ ray_tet_patch_intersect, \
+ tet_patch_centroid, \
+ tet_patch_bbox
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
@@ -131,10 +131,10 @@
self._set_up_patches(vertices, indices)
elif self.num_verts_per_elem == 10:
self.primitives = malloc(self.num_prim * sizeof(TetPatch))
- self.get_centroid = tetPatch_centroid
- self.get_bbox = tetPatch_bbox
- self.get_intersect = ray_tetPatch_intersect
- self._set_up_tetPatches(vertices, indices)
+ self.get_centroid = tet_patch_centroid
+ self.get_bbox = tet_patch_bbox
+ self.get_intersect = ray_tet_patch_intersect
+ self._set_up_tet_patches(vertices, indices)
else:
self.primitives = malloc(self.num_prim * sizeof(Triangle))
self.get_centroid = triangle_centroid
@@ -173,22 +173,22 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void _set_up_tetPatches(self, np.float64_t[:, :] vertices,
+ cdef void _set_up_tet_patches(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:
- cdef TetPatch* tetPatch
+ cdef TetPatch* tet_patch
cdef np.int64_t i, j, k, ind, idim
cdef np.int64_t offset, prim_index
for i in range(self.num_elem):
offset = self.num_prim_per_elem*i
for j in range(self.num_prim_per_elem): # for each face
prim_index = offset + j
- tetPatch = &( <TetPatch*> self.primitives)[prim_index]
+ tet_patch = &( <TetPatch*> self.primitives)[prim_index]
self.prim_ids[prim_index] = prim_index
- tetPatch.elem_id = i
+ tet_patch.elem_id = i
for k in range(6): # for each vertex
ind = tet10_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
- tetPatch.v[k][idim] = vertices[indices[i, ind]][idim]
+ tet_patch.v[k][idim] = vertices[indices[i, ind]][idim]
self.get_centroid(self.primitives,
prim_index,
self.centroids[prim_index])
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/lib/mesh_triangulation.h
--- a/yt/utilities/lib/mesh_triangulation.h
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -67,8 +67,8 @@
// This is used to select faces from a second-order tet element
int tet10_faces[4][6] = {
- {0, 1, 3, 7, 4, 8},
- {2, 3, 1, 5, 9, 8},
- {0, 3, 2, 6, 7, 9},
- {0, 2, 1, 4, 6, 5}
+ {0, 1, 3, 4, 8, 7},
+ {2, 3, 1, 9, 8, 5},
+ {0, 3, 2, 7, 9, 6},
+ {0, 2, 1, 6, 5, 4}
};
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/lib/primitives.pxd
--- a/yt/utilities/lib/primitives.pxd
+++ b/yt/utilities/lib/primitives.pxd
@@ -83,33 +83,33 @@
np.float64_t[6][3] v # 6 vertices per patch
np.int64_t elem_id
-cdef RayHitData compute_tetPatch_hit(cython.floating[6][3] verts,
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil
-cdef void tetPatchSurfaceFunc(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] S) nogil
-cdef void tetPatchSurfaceDerivU(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] Su) nogil
-cdef void tetPatchSurfaceDerivV(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil
-cdef np.int64_t ray_tetPatch_intersect(const void* primitives,
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil
-cdef void tetPatch_centroid(const void *primitives,
+cdef void tet_patch_centroid(const void *primitives,
const np.int64_t item,
np.float64_t[3] centroid) nogil
-cdef void tetPatch_bbox(const void *primitives,
+cdef void tet_patch_bbox(const void *primitives,
const np.int64_t item,
BBox* bbox) nogil
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -305,7 +305,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void tetPatchSurfaceFunc(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] S) nogil:
@@ -323,7 +323,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void tetPatchSurfaceDerivU(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] Su) nogil:
@@ -340,7 +340,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void tetPatchSurfaceDerivV(const cython.floating[6][3] verts,
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil:
@@ -358,7 +358,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef RayHitData compute_tetPatch_hit(cython.floating[6][3] verts,
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil:
@@ -385,7 +385,7 @@
cdef cython.floating u = 0.0
cdef cython.floating v = 0.0
cdef cython.floating[3] S
- tetPatchSurfaceFunc(verts, u, v, S)
+ tet_patchSurfaceFunc(verts, u, v, S)
cdef cython.floating fu = dot(N1, S) + d1
cdef cython.floating fv = dot(N2, S) + d2
cdef cython.floating err = fmax(fabs(fu), fabs(fv))
@@ -399,8 +399,8 @@
cdef cython.floating J11, J12, J21, J22, det
while ((err > tol) and (iterations < max_iter)):
# compute the Jacobian
- tetPatchSurfaceDerivU(verts, u, v, Su)
- tetPatchSurfaceDerivV(verts, u, v, Sv)
+ tet_patchSurfaceDerivU(verts, u, v, Su)
+ tet_patchSurfaceDerivV(verts, u, v, Sv)
J11 = dot(N1, Su)
J12 = dot(N1, Sv)
J21 = dot(N2, Su)
@@ -411,7 +411,7 @@
u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
- tetPatchSurfaceFunc(verts, u, v, S)
+ tet_patchSurfaceFunc(verts, u, v, S)
fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
@@ -433,13 +433,13 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef np.int64_t ray_tetPatch_intersect(const void* primitives,
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil:
- cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
- cdef RayHitData hd = compute_tetPatch_hit(tetPatch.v, ray.origin, ray.direction)
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.origin, ray.direction)
# only count this is it's the closest hit
if (hd.t < ray.t_near or hd.t > ray.t_far):
@@ -448,7 +448,7 @@
if ((0 <= hd.u <= 1.0) and (0 <= hd.v <= 1.0) and hd.converged):
# we have a hit, so update ray information
ray.t_far = hd.t
- ray.elem_id = tetPatch.elem_id
+ ray.elem_id = tet_patch.elem_id
return True
return False
@@ -457,19 +457,19 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void tetPatch_centroid(const void *primitives,
+cdef void tet_patch_centroid(const void *primitives,
const np.int64_t item,
np.float64_t[3] centroid) nogil:
cdef np.int64_t i, j
- cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
for j in range(3):
centroid[j] = 0.0
for i in range(6):
for j in range(3):
- centroid[j] += tetPatch.v[i][j]
+ centroid[j] += tet_patch.v[i][j]
for j in range(3):
centroid[j] /= 6.0
@@ -478,18 +478,18 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void tetPatch_bbox(const void *primitives,
+cdef void tet_patch_bbox(const void *primitives,
const np.int64_t item,
BBox* bbox) nogil:
cdef np.int64_t i, j
- cdef TetPatch tetPatch = (<TetPatch*> primitives)[item]
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
for j in range(3):
- bbox.left_edge[j] = tetPatch.v[0][j]
- bbox.right_edge[j] = tetPatch.v[0][j]
+ bbox.left_edge[j] = tet_patch.v[0][j]
+ bbox.right_edge[j] = tet_patch.v[0][j]
for i in range(1, 6):
for j in range(3):
- bbox.left_edge[j] = fmin(bbox.left_edge[j], tetPatch.v[i][j])
- bbox.right_edge[j] = fmax(bbox.right_edge[j], tetPatch.v[i][j])
+ bbox.left_edge[j] = fmin(bbox.left_edge[j], tet_patch.v[i][j])
+ bbox.right_edge[j] = fmax(bbox.right_edge[j], tet_patch.v[i][j])
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/utilities/orientation.py
--- a/yt/utilities/orientation.py
+++ b/yt/utilities/orientation.py
@@ -23,7 +23,7 @@
def _aligned(a, b):
aligned_component = np.abs(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
return np.isclose(aligned_component, 1.0, 1.0e-13)
-
+
def _validate_unit_vectors(normal_vector, north_vector):
@@ -35,7 +35,6 @@
if not np.dot(normal_vector, normal_vector) > 0:
raise YTException("normal_vector cannot be the zero vector.")
-
if north_vector is not None and _aligned(north_vector, normal_vector):
raise YTException("normal_vector and north_vector cannot be aligned.")
@@ -85,7 +84,7 @@
t = np.cross(normal_vector, vecs).sum(axis=1)
ax = t.argmax()
east_vector = np.cross(vecs[ax, :], normal_vector).ravel()
- # self.north_vector must remain None otherwise rotations about a fixed axis will break.
+ # self.north_vector must remain None otherwise rotations about a fixed axis will break.
# The north_vector calculated here will still be included in self.unit_vectors.
north_vector = np.cross(normal_vector, east_vector).ravel()
else:
diff -r 05262ce494229229ad11ec309695f445c08f0b8c -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -158,7 +158,7 @@
def position():
doc = '''
- The location of the camera.
+ The location of the camera.
Parameters
----------
@@ -483,11 +483,11 @@
>>> sc = Scene()
>>> cam = sc.add_camera()
>>> # rotate the camera by pi / 4 radians:
- >>> cam.rotate(np.pi/4.0)
+ >>> cam.rotate(np.pi/4.0)
>>> # rotate the camera about the y-axis instead of cam.north_vector:
- >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
>>> # rotate the camera about the origin instead of its own position:
- >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
"""
rotate_all = rot_vector is not None
@@ -593,7 +593,7 @@
>>> sc = Scene()
>>> cam = sc.add_camera(ds)
>>> # roll the camera by pi / 4 radians:
- >>> cam.roll(np.pi/4.0)
+ >>> cam.roll(np.pi/4.0)
>>> # roll the camera about the origin instead of its own position:
>>> cam.roll(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
@@ -627,7 +627,7 @@
>>> import yt
>>> import numpy as np
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
>>> cam = sc.camera
>>> for i in cam.iter_rotate(np.pi, 10):
@@ -688,7 +688,7 @@
def zoom(self, factor):
r"""Change the width of the FOV of the camera.
- This will appear to zoom the camera in by some `factor` toward the
+ This will appear to zoom the camera in by some `factor` toward the
focal point along the current view direction, but really it's just
changing the width of the field of view.
https://bitbucket.org/yt_analysis/yt/commits/cec7d2011495/
Changeset: cec7d2011495
Branch: yt
User: al007
Date: 2016-09-30 21:30:17+00:00
Summary: Merge yt-tip.
Affected #: 61 files
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 doc/source/analyzing/generating_processed_data.rst
--- a/doc/source/analyzing/generating_processed_data.rst
+++ b/doc/source/analyzing/generating_processed_data.rst
@@ -26,13 +26,16 @@
sizes into a fixed-size array that appears like an image. This process is that
of pixelization, which yt handles transparently internally. You can access
this functionality by constructing a
-:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` (or
-:class:`~yt.visualization.fixed_resolution.ObliqueFixedResolutionBuffer`) and
-supplying to it your :class:`~yt.data_objects.data_containers.YTSelectionContainer2D`
+:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` and supplying
+to it your :class:`~yt.data_objects.data_containers.YTSelectionContainer2D`
object, as well as some information about how you want the final image to look.
You can specify both the bounds of the image (in the appropriate x-y plane) and
-the resolution of the output image. You can then have yt pixelize any
-field you like.
+the resolution of the output image. You can then have yt pixelize any field
+you like.
+
+.. note:: In previous versions of yt, there was a special class of
+ FixedResolutionBuffer for off-axis slices. This is no longer
+ necessary.
To create :class:`~yt.data_objects.data_containers.YTSelectionContainer2D` objects, you can
access them as described in :ref:`data-objects`, specifically the section
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 doc/source/analyzing/objects.rst
--- a/doc/source/analyzing/objects.rst
+++ b/doc/source/analyzing/objects.rst
@@ -312,9 +312,12 @@
| easily lead to empty data for non-intersecting regions.
| Usage: ``slice(axis, coord, ds, data_source=sph)``
-**Boolean Regions**
- | **Note: not yet implemented in yt 3.0**
- | Usage: ``boolean()``
+**Union Regions**
+ | Usage: ``union()``
+ | See :ref:`boolean_data_objects`.
+
+**Intersection Regions**
+ | Usage: ``intersection()``
| See :ref:`boolean_data_objects`.
**Filter**
@@ -614,37 +617,48 @@
Combining Objects: Boolean Data Objects
---------------------------------------
-.. note:: Boolean Data Objects have not yet been ported to yt 3.0 from
- yt 2.x. If you are interested in aiding in this port, please contact
- the yt-dev mailing list. Until it is ported, this functionality below
- will not work.
+A special type of data object is the *boolean* data object, which works with
+three-dimensional data selection. It is built by relating already existing
+data objects with the bitwise operators for AND, OR and XOR, as well as the
+subtraction operator. These are created by using the operators ``&`` for an
+intersection ("AND"), ``|`` for a union ("OR"), ``^`` for an exclusive or
+("XOR"), and ``+`` and ``-`` for addition ("OR") and subtraction ("NEG").
+Here are some examples:
-A special type of data object is the *boolean* data object.
-It works only on three-dimensional objects.
-It is built by relating already existing data objects with boolean operators.
-The boolean logic may be nested using parentheses, and
-it supports the standard "AND", "OR", and "NOT" operators:
+.. code-block:: python
-* **"AND"** When two data objects are related with an "AND", the combined
- data object is the volume of the simulation covered by both objects, and
- not by just a single object.
-* **"OR"** When two data objects are related with an "OR", the combined
- data object is the volume(s) of the simulation covered by either of the
- objects.
- For example, this may be used to combine disjoint objects into one.
-* **"NOT"** When two data objects are related with a "NOT", the combined
- data object is the volume of the first object that the second does not
- cover.
- For example, this may be used to cut out part(s) of the first data object
- utilizing the second data object.
-* **"(" or ")"** Nested logic is surrounded by parentheses. The order of
- operations is such that the boolean logic is evaluated inside the
- inner-most parentheses, first, then goes upwards.
- The logic is read left-to-right at all levels (crucial for the "NOT"
- operator).
+ import yt
+ ds = yt.load("snapshot_010.hdf5")
-Please see the :ref:`cookbook` for some examples of how to use the boolean
-data object.
+ sp1 = ds.sphere("c", (0.1, "unitary"))
+ sp2 = ds.sphere(sp1.center + 2.0 * sp1.radius, (0.2, "unitary"))
+ sp3 = ds.sphere("c", (0.05, "unitary"))
+
+ new_obj = sp1 + sp2
+ cutout = sp1 - sp3
+ sp4 = sp1 ^ sp2
+ sp5 = sp1 & sp2
+
+
+Note that the ``+`` operation and the ``|`` operation are identical. For when
+multiple objects are to be combined in an intersection or a union, there are
+the data objects ``intersection`` and ``union`` which can be called, and which
+will yield slightly higher performance than a sequence of calls to ``+`` or
+``&``. For instance:
+
+.. code-block:: python
+
+ import yt
+ ds = yt.load("Enzo_64/DD0043/data0043")
+ sp1 = ds.sphere( (0.1, 0.2, 0.3), (0.05, "unitary"))
+ sp2 = ds.sphere( (0.2, 0.2, 0.3), (0.10, "unitary"))
+ sp3 = ds.sphere( (0.3, 0.2, 0.3), (0.15, "unitary"))
+
+ isp = ds.intersection( [sp1, sp2, sp3] )
+ usp = ds.union( [sp1, sp2, sp3] )
+
+The ``isp`` and ``usp`` objects will act the same as a set of chained ``&`` and
+``|`` operations (respectively) but are somewhat easier to construct.
.. _extracting-connected-sets:
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -1264,23 +1264,37 @@
.. code-block:: python
- import yt
- import numpy
- from yt.utilities.exodusII_reader import get_data
+ import yt
+ import numpy as np
- coords, connectivity, data = get_data("MOOSE_sample_data/out.e-s010")
+ coords = np.array([[0.0, 0.0],
+ [1.0, 0.0],
+ [1.0, 1.0],
+ [0.0, 1.0]], dtype=np.float64)
-This uses a publically available `MOOSE <http://mooseframework.org/>`
-dataset along with the get_data function to parse the coords, connectivity,
-and data. Then, these can be loaded as an in-memory dataset as follows:
+ connect = np.array([[0, 1, 3],
+ [1, 2, 3]], dtype=np.int64)
+
+ data = {}
+ data['connect1', 'test'] = np.array([[0.0, 1.0, 3.0],
+ [1.0, 2.0, 3.0]], dtype=np.float64)
+
+Here, we have made up a simple, 2D unstructured mesh dataset consisting of two
+triangles and one node-centered data field. This data can be loaded as an in-memory
+dataset as follows:
.. code-block:: python
- mesh_id = 0
- ds = yt.load_unstructured_mesh(data[mesh_id], connectivity[mesh_id], coords[mesh_id])
+ ds = yt.load_unstructured_mesh(connect, coords, data)
-Note that load_unstructured_mesh can take either a single or a list of meshes.
-Here, we have selected only the first mesh to load.
+Note that load_unstructured_mesh can take either a single mesh or a list of meshes.
+Here, we only have one mesh. The in-memory dataset can then be visualized as usual,
+e.g.:
+
+.. code-block:: python
+
+ sl = yt.SlicePlot(ds, 'z', 'test')
+ sl.annotate_mesh_lines()
.. rubric:: Caveats
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 doc/source/reference/api/api.rst
--- a/doc/source/reference/api/api.rst
+++ b/doc/source/reference/api/api.rst
@@ -49,7 +49,6 @@
~yt.visualization.fixed_resolution.FixedResolutionBuffer
~yt.visualization.fixed_resolution.ParticleImageBuffer
~yt.visualization.fixed_resolution.CylindricalFixedResolutionBuffer
- ~yt.visualization.fixed_resolution.ObliqueFixedResolutionBuffer
~yt.visualization.fixed_resolution.OffAxisProjectionFixedResolutionBuffer
Data Sources
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -39,7 +39,7 @@
local_owls_000:
- yt/frontends/owls/tests/test_outputs.py
- local_pw_006:
+ local_pw_008:
- yt/visualization/tests/test_plotwindow.py:test_attributes
- yt/visualization/tests/test_plotwindow.py:test_attributes_wt
- yt/visualization/tests/test_profile_plots.py:test_phase_plot_attributes
@@ -50,7 +50,7 @@
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_004:
+ local_varia_005:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -49,7 +49,9 @@
YTFieldNotFound, \
YTFieldTypeNotFound, \
YTDataSelectorNotImplemented, \
- YTDimensionalityError
+ YTDimensionalityError, \
+ YTBooleanObjectError, \
+ YTBooleanObjectsWrongDataset
from yt.utilities.lib.marching_cubes import \
march_cubes_grid, march_cubes_grid_flux
from yt.utilities.parallel_tools.parallel_analysis_interface import \
@@ -1415,7 +1417,6 @@
self.index._identify_base_chunk(self)
return self._current_chunk.fcoords_vertex
-
class YTSelectionContainer0D(YTSelectionContainer):
_spatial = False
_dimensionality = 0
@@ -1872,6 +1873,86 @@
"""
return self.quantities.total_quantity(("index", "cell_volume"))
+ def __or__(self, other):
+ if not isinstance(other, YTSelectionContainer3D):
+ raise YTBooleanObjectError(other)
+ if self.ds is not other.ds:
+ raise YTBooleanObjectsWrongDataset()
+ # Should maybe do something with field parameters here
+ return YTBooleanContainer("OR", self, other, ds = self.ds)
+
+ def __invert__(self):
+ # ~obj
+ asel = yt.geometry.selection_routines.AlwaysSelector(self.ds)
+ return YTBooleanContainer("NOT", self, asel, ds = self.ds)
+
+ def __xor__(self, other):
+ if not isinstance(other, YTSelectionContainer3D):
+ raise YTBooleanObjectError(other)
+ if self.ds is not other.ds:
+ raise YTBooleanObjectsWrongDataset()
+ return YTBooleanContainer("XOR", self, other, ds = self.ds)
+
+ def __and__(self, other):
+ if not isinstance(other, YTSelectionContainer3D):
+ raise YTBooleanObjectError(other)
+ if self.ds is not other.ds:
+ raise YTBooleanObjectsWrongDataset()
+ return YTBooleanContainer("AND", self, other, ds = self.ds)
+
+ def __add__(self, other):
+ return self.__or__(other)
+
+ def __sub__(self, other):
+ if not isinstance(other, YTSelectionContainer3D):
+ raise YTBooleanObjectError(other)
+ if self.ds is not other.ds:
+ raise YTBooleanObjectsWrongDataset()
+ return YTBooleanContainer("NEG", self, other, ds = self.ds)
+
+class YTBooleanContainer(YTSelectionContainer3D):
+ """
+ This is a boolean operation, accepting AND, OR, XOR, and NOT for combining
+ multiple data objects.
+
+ This object is not designed to be created directly; it is designed to be
+ created implicitly by using one of the bitwise operations (&, |, ^, ~) on
+ one or two other data objects. These correspond to the appropriate boolean
+ operations, and the resultant object can be nested.
+
+ Parameters
+ ----------
+ op : string
+ Can be AND, OR, XOR, NOT or NEG.
+ dobj1 : YTSelectionContainer3D
+ The first selection object
+ dobj2 : YTSelectionContainer3D
+ The second object
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+ >>> sp = ds.sphere("c", 0.1)
+ >>> dd = ds.r[:,:,:]
+ >>> new_obj = sp ^ dd
+ >>> print(new_obj.sum("cell_volume"), dd.sum("cell_volume") -
+ ... sp.sum("cell_volume"))
+ """
+ _type_name = "bool"
+ _con_args = ("op", "dobj1", "dobj2")
+ def __init__(self, op, dobj1, dobj2, ds = None, field_parameters = None,
+ data_source = None):
+ YTSelectionContainer3D.__init__(self, None, ds, field_parameters,
+ data_source)
+ self.op = op.upper()
+ self.dobj1 = dobj1
+ self.dobj2 = dobj2
+ name = "Boolean%sSelector" % (self.op,)
+ sel_cls = getattr(yt.geometry.selection_routines, name)
+ self._selector = sel_cls(self)
+
# Many of these items are set up specifically to ensure that
# we are not breaking old pickle files. This means we must only call the
# _reconstruct_object and that we cannot mandate any additional arguments to
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -463,12 +463,12 @@
self.fields = ensure_list(fields) + [k for k in self.field_data.keys()
if k not in self._key_fields]
from yt.visualization.plot_window import get_oblique_window_parameters, PWViewerMPL
- from yt.visualization.fixed_resolution import ObliqueFixedResolutionBuffer
+ from yt.visualization.fixed_resolution import FixedResolutionBuffer
(bounds, center_rot) = get_oblique_window_parameters(normal, center, width, self.ds)
pw = PWViewerMPL(
self, bounds, fields=self.fields, origin='center-window',
periodic=False, oblique=True,
- frb_generator=ObliqueFixedResolutionBuffer,
+ frb_generator=FixedResolutionBuffer,
plot_type='OffAxisSlice')
if axes_unit is not None:
pw.set_axes_unit(axes_unit)
@@ -476,8 +476,8 @@
return pw
def to_frb(self, width, resolution, height=None, periodic=False):
- r"""This function returns an ObliqueFixedResolutionBuffer generated
- from this object.
+ r"""This function returns a FixedResolutionBuffer generated from this
+ object.
An ObliqueFixedResolutionBuffer is an object that accepts a
variable-resolution 2D object and transforms it into an NxM bitmap that
@@ -526,9 +526,9 @@
height = self.ds.quan(height[0], height[1])
if not iterable(resolution):
resolution = (resolution, resolution)
- from yt.visualization.fixed_resolution import ObliqueFixedResolutionBuffer
+ from yt.visualization.fixed_resolution import FixedResolutionBuffer
bounds = (-width/2.0, width/2.0, -height/2.0, height/2.0)
- frb = ObliqueFixedResolutionBuffer(self, bounds, resolution,
+ frb = FixedResolutionBuffer(self, bounds, resolution,
periodic=periodic)
return frb
@@ -869,3 +869,75 @@
@property
def fwidth(self):
return self.base_object.fwidth[self._cond_ind,:]
+
+class YTIntersectionContainer3D(YTSelectionContainer3D):
+ """
+ This is a more efficient method of selecting the intersection of multiple
+ data selection objects.
+
+ Creating one of these objects returns the intersection of all of the
+ sub-objects; it is designed to be a faster method than chaining & ("and")
+ operations to create a single, large intersection.
+
+ Parameters
+ ----------
+ data_objects : Iterable of YTSelectionContainer3D
+ The data objects to intersect
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load("RedshiftOutput0005")
+ >>> sp1 = ds.sphere((0.4, 0.5, 0.6), 0.15)
+ >>> sp2 = ds.sphere((0.38, 0.51, 0.55), 0.1)
+ >>> sp3 = ds.sphere((0.35, 0.5, 0.6), 0.15)
+ >>> new_obj = ds.intersection((sp1, sp2, sp3))
+ >>> print(new_obj.sum("cell_volume"))
+ """
+ _type_name = "intersection"
+ _con_args = ("data_objects",)
+ def __init__(self, data_objects, ds = None, field_parameters = None,
+ data_source = None):
+ YTSelectionContainer3D.__init__(self, None, ds, field_parameters,
+ data_source)
+ # ensure_list doesn't check for tuples
+ if isinstance(data_objects, tuple):
+ data_objects = list(data_objects)
+ self.data_objects = ensure_list(data_objects)
+
+class YTDataObjectUnion(YTSelectionContainer3D):
+ """
+ This is a more efficient method of selecting the union of multiple
+ data selection objects.
+
+ Creating one of these objects returns the union of all of the sub-objects;
+ it is designed to be a faster method than chaining | (or) operations to
+ create a single, large union.
+
+ Parameters
+ ----------
+ data_objects : Iterable of YTSelectionContainer3D
+ The data objects to union
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+ >>> sp1 = ds.sphere((0.4, 0.5, 0.6), 0.1)
+ >>> sp2 = ds.sphere((0.3, 0.5, 0.15), 0.1)
+ >>> sp3 = ds.sphere((0.5, 0.5, 0.9), 0.1)
+ >>> new_obj = ds.union((sp1, sp2, sp3))
+ >>> print(new_obj.sum("cell_volume"))
+ """
+ _type_name = "union"
+ _con_args = ("data_objects",)
+ def __init__(self, data_objects, ds = None, field_parameters = None,
+ data_source = None):
+ YTSelectionContainer3D.__init__(self, None, ds, field_parameters,
+ data_source)
+ # ensure_list doesn't check for tuples
+ if isinstance(data_objects, tuple):
+ data_objects = list(data_objects)
+ self.data_objects = ensure_list(data_objects)
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -20,6 +20,7 @@
import os
import time
import weakref
+import warnings
from collections import defaultdict
from yt.extern.six import add_metaclass, string_types
@@ -214,7 +215,7 @@
obj = _cached_datasets[apath]
return obj
- def __init__(self, filename, dataset_type=None, file_style=None,
+ def __init__(self, filename, dataset_type=None, file_style=None,
units_override=None, unit_system="cgs"):
"""
Base class for generating new output types. Principally consists of
@@ -338,7 +339,7 @@
in that directory, and a list of subdirectories. It should return a
list of filenames (defined relative to the supplied directory) and a
boolean as to whether or not further directories should be recursed.
-
+
This function doesn't need to catch all possibilities, nor does it need
to filter possibilities.
"""
@@ -938,7 +939,7 @@
"dangerous option that may yield inconsistent results, and must be "
"used very carefully, and only if you know what you want from it.")
for unit, cgs in [("length", "cm"), ("time", "s"), ("mass", "g"),
- ("velocity","cm/s"), ("magnetic","gauss"),
+ ("velocity","cm/s"), ("magnetic","gauss"),
("temperature","K")]:
val = self.units_override.get("%s_unit" % unit, None)
if val is not None:
@@ -1043,7 +1044,7 @@
self._quan = functools.partial(YTQuantity, registry=self.unit_registry)
return self._quan
- def add_field(self, name, function=None, **kwargs):
+ def add_field(self, name, function=None, sampling_type=None, **kwargs):
"""
Dataset-specific call to add_field
@@ -1081,7 +1082,18 @@
if not override and name in self.field_info:
mylog.warning("Field %s already exists. To override use " +
"force_override=True.", name)
- self.field_info.add_field(name, function=function, **kwargs)
+ if kwargs.setdefault('particle_type', False):
+ if sampling_type is not None and sampling_type != "particle":
+ raise RuntimeError("Clashing definition of 'sampling_type' and "
+ "'particle_type'. Note that 'particle_type' is "
+ "deprecated. Please just use 'sampling_type'.")
+ else:
+ sampling_type = "particle"
+ if sampling_type is None:
+ warnings.warn("Because 'sampling_type' not specified, yt will "
+ "assume a cell 'sampling_type'")
+ sampling_type = "cell"
+ self.field_info.add_field(name, sampling_type, function=function, **kwargs)
self.field_info._show_field_errors.append(name)
deps, _ = self.field_info.check_derived_fields([name])
self.field_dependencies.update(deps)
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/data_objects/tests/test_boolean_regions.py
--- /dev/null
+++ b/yt/data_objects/tests/test_boolean_regions.py
@@ -0,0 +1,429 @@
+from yt.testing import \
+ fake_amr_ds, \
+ assert_array_equal
+import numpy as np
+
+# We use morton indices in this test because they are single floating point
+# values that uniquely identify each cell. That's a convenient way to compare
+# inclusion in set operations, since there are no duplicates.
+
+def test_boolean_spheres_no_overlap():
+ r"""Test to make sure that boolean objects (spheres, no overlap)
+ behave the way we expect.
+
+ Test non-overlapping spheres. This also checks that the original spheres
+ don't change as part of constructing the booleans.
+ """
+ ds = fake_amr_ds()
+ sp1 = ds.sphere([0.25, 0.25, 0.25], 0.15)
+ sp2 = ds.sphere([0.75, 0.75, 0.75], 0.15)
+ # Store the original indices
+ i1 = sp1["index","morton_index"]
+ i1.sort()
+ i2 = sp2["index","morton_index"]
+ i2.sort()
+ ii = np.concatenate((i1, i2))
+ ii.sort()
+ # Make some booleans
+ bo1 = sp1 & sp2
+ bo2 = sp1 - sp2
+ bo3 = sp1 | sp2 # also works with +
+ bo4 = ds.union([sp1, sp2])
+ bo5 = ds.intersection([sp1, sp2])
+ # This makes sure the original containers didn't change.
+ new_i1 = sp1["index","morton_index"]
+ new_i1.sort()
+ new_i2 = sp2["index","morton_index"]
+ new_i2.sort()
+ assert_array_equal(new_i1, i1)
+ assert_array_equal(new_i2, i2)
+ # Now make sure the indices also behave as we expect.
+ empty = np.array([])
+ assert_array_equal(bo1["index","morton_index"], empty)
+ assert_array_equal(bo5["index","morton_index"], empty)
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ assert_array_equal(b2, i1)
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b3, ii)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ assert_array_equal(b4, ii)
+ bo6 = sp1 ^ sp2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_spheres_overlap():
+ r"""Test to make sure that boolean objects (spheres, overlap)
+ behave the way we expect.
+
+ Test overlapping spheres.
+ """
+ ds = fake_amr_ds()
+ sp1 = ds.sphere([0.45, 0.45, 0.45], 0.15)
+ sp2 = ds.sphere([0.55, 0.55, 0.55], 0.15)
+ # Get indices of both.
+ i1 = sp1["index","morton_index"]
+ i2 = sp2["index","morton_index"]
+ # Make some booleans
+ bo1 = sp1 & sp2
+ bo2 = sp1 - sp2
+ bo3 = sp1 | sp2
+ bo4 = ds.union([sp1, sp2])
+ bo5 = ds.intersection([sp1, sp2])
+ # Now make sure the indices also behave as we expect.
+ lens = np.intersect1d(i1, i2)
+ apple = np.setdiff1d(i1, i2)
+ both = np.union1d(i1, i2)
+ b1 = bo1["index","morton_index"]
+ b1.sort()
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b1, lens)
+ assert_array_equal(b2, apple)
+ assert_array_equal(b3, both)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ assert_array_equal(b1, b5)
+ bo6 = sp1 ^ sp2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_regions_no_overlap():
+ r"""Test to make sure that boolean objects (regions, no overlap)
+ behave the way we expect.
+
+ Test non-overlapping regions. This also checks that the original regions
+ don't change as part of constructing the booleans.
+ """
+ ds = fake_amr_ds()
+ re1 = ds.region([0.25]*3, [0.2]*3, [0.3]*3)
+ re2 = ds.region([0.65]*3, [0.6]*3, [0.7]*3)
+ # Store the original indices
+ i1 = re1["index","morton_index"]
+ i1.sort()
+ i2 = re2["index","morton_index"]
+ i2.sort()
+ ii = np.concatenate((i1, i2))
+ ii.sort()
+ # Make some booleans
+ bo1 = re1 & re2
+ bo2 = re1 - re2
+ bo3 = re1 | re2
+ bo4 = ds.union([re1, re2])
+ bo5 = ds.intersection([re1, re2])
+ # This makes sure the original containers didn't change.
+ new_i1 = re1["index","morton_index"]
+ new_i1.sort()
+ new_i2 = re2["index","morton_index"]
+ new_i2.sort()
+ assert_array_equal(new_i1, i1)
+ assert_array_equal(new_i2, i2)
+ # Now make sure the indices also behave as we expect.
+ empty = np.array([])
+ assert_array_equal(bo1["index","morton_index"], empty)
+ assert_array_equal(bo5["index","morton_index"], empty)
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ assert_array_equal(b2, i1 )
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b3, ii)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ bo6 = re1 ^ re2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_regions_overlap():
+ r"""Test to make sure that boolean objects (regions, overlap)
+ behave the way we expect.
+
+ Test overlapping regions.
+ """
+ ds = fake_amr_ds()
+ re1 = ds.region([0.55]*3, [0.5]*3, [0.6]*3)
+ re2 = ds.region([0.6]*3, [0.55]*3, [0.65]*3)
+ # Get indices of both.
+ i1 = re1["index","morton_index"]
+ i2 = re2["index","morton_index"]
+ # Make some booleans
+ bo1 = re1 & re2
+ bo2 = re1 - re2
+ bo3 = re1 | re2
+ bo4 = ds.union([re1, re2])
+ bo5 = ds.intersection([re1, re2])
+ # Now make sure the indices also behave as we expect.
+ cube = np.intersect1d(i1, i2)
+ bite_cube = np.setdiff1d(i1, i2)
+ both = np.union1d(i1, i2)
+ b1 = bo1["index","morton_index"]
+ b1.sort()
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b1, cube)
+ assert_array_equal(b2, bite_cube)
+ assert_array_equal(b3, both)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ assert_array_equal(b1, b5)
+ bo6 = re1 ^ re2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_cylinders_no_overlap():
+ r"""Test to make sure that boolean objects (cylinders, no overlap)
+ behave the way we expect.
+
+ Test non-overlapping cylinders. This also checks that the original cylinders
+ don't change as part of constructing the booleans.
+ """
+ ds = fake_amr_ds()
+ cyl1 = ds.disk([0.25]*3, [1, 0, 0], 0.1, 0.1)
+ cyl2 = ds.disk([0.75]*3, [1, 0, 0], 0.1, 0.1)
+ # Store the original indices
+ i1 = cyl1["index","morton_index"]
+ i1.sort()
+ i2 = cyl2["index","morton_index"]
+ i2.sort()
+ ii = np.concatenate((i1, i2))
+ ii.sort()
+ # Make some booleans
+ bo1 = cyl1 & cyl2
+ bo2 = cyl1 - cyl2
+ bo3 = cyl1 | cyl2
+ bo4 = ds.union([cyl1, cyl2])
+ bo5 = ds.intersection([cyl1, cyl2])
+ # This makes sure the original containers didn't change.
+ new_i1 = cyl1["index","morton_index"]
+ new_i1.sort()
+ new_i2 = cyl2["index","morton_index"]
+ new_i2.sort()
+ assert_array_equal(new_i1, i1)
+ assert_array_equal(new_i2, i2)
+ # Now make sure the indices also behave as we expect.
+ empty = np.array([])
+ assert_array_equal(bo1["index","morton_index"], empty)
+ assert_array_equal(bo5["index","morton_index"], empty)
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ assert_array_equal(b2, i1)
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b3, ii)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ bo6 = cyl1 ^ cyl2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_cylinders_overlap():
+ r"""Test to make sure that boolean objects (cylinders, overlap)
+ behave the way we expect.
+
+ Test overlapping cylinders.
+ """
+ ds = fake_amr_ds()
+ cyl1 = ds.disk([0.45]*3, [1, 0, 0], 0.2, 0.2)
+ cyl2 = ds.disk([0.55]*3, [1, 0, 0], 0.2, 0.2)
+ # Get indices of both.
+ i1 = cyl1["index","morton_index"]
+ i2 = cyl2["index","morton_index"]
+ # Make some booleans
+ bo1 = cyl1 & cyl2
+ bo2 = cyl1 - cyl2
+ bo3 = cyl1 | cyl2
+ bo4 = ds.union([cyl1, cyl2])
+ bo5 = ds.intersection([cyl1, cyl2])
+ # Now make sure the indices also behave as we expect.
+ vlens = np.intersect1d(i1, i2)
+ bite_disk = np.setdiff1d(i1, i2)
+ both = np.union1d(i1, i2)
+ b1 = bo1["index","morton_index"]
+ b1.sort()
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b1, vlens)
+ assert_array_equal(b2, bite_disk)
+ assert_array_equal(b3, both)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ assert_array_equal(b1, b5)
+ bo6 = cyl1 ^ cyl2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_ellipsoids_no_overlap():
+ r"""Test to make sure that boolean objects (ellipsoids, no overlap)
+ behave the way we expect.
+
+ Test non-overlapping ellipsoids. This also checks that the original
+ ellipsoids don't change as part of constructing the booleans.
+ """
+ ds = fake_amr_ds()
+ ell1 = ds.ellipsoid([0.25]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ ell2 = ds.ellipsoid([0.75]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ # Store the original indices
+ i1 = ell1["index","morton_index"]
+ i1.sort()
+ i2 = ell2["index","morton_index"]
+ i2.sort()
+ ii = np.concatenate((i1, i2))
+ ii.sort()
+ # Make some booleans
+ bo1 = ell1 & ell2
+ bo2 = ell1 - ell2
+ bo3 = ell1 | ell2
+ bo4 = ds.union([ell1, ell2])
+ bo5 = ds.intersection([ell1, ell2])
+ # This makes sure the original containers didn't change.
+ new_i1 = ell1["index","morton_index"]
+ new_i1.sort()
+ new_i2 = ell2["index","morton_index"]
+ new_i2.sort()
+ assert_array_equal(new_i1, i1 )
+ assert_array_equal(new_i2, i2)
+ # Now make sure the indices also behave as we expect.
+ empty = np.array([])
+ assert_array_equal(bo1["index","morton_index"], empty)
+ assert_array_equal(bo5["index","morton_index"], empty)
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ assert_array_equal(b2, i1)
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b3, ii)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ bo6 = ell1 ^ ell2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_ellipsoids_overlap():
+ r"""Test to make sure that boolean objects (ellipsoids, overlap)
+ behave the way we expect.
+
+ Test overlapping ellipsoids.
+ """
+ ds = fake_amr_ds()
+ ell1 = ds.ellipsoid([0.45]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ ell2 = ds.ellipsoid([0.55]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ # Get indices of both.
+ i1 = ell1["index","morton_index"]
+ i2 = ell2["index","morton_index"]
+ # Make some booleans
+ bo1 = ell1 & ell2
+ bo2 = ell1 - ell2
+ bo3 = ell1 | ell2
+ bo4 = ds.union([ell1, ell2])
+ bo5 = ds.intersection([ell1, ell2])
+ # Now make sure the indices also behave as we expect.
+ overlap = np.intersect1d(i1, i2)
+ diff = np.setdiff1d(i1, i2)
+ both = np.union1d(i1, i2)
+ b1 = bo1["index","morton_index"]
+ b1.sort()
+ b2 = bo2["index","morton_index"]
+ b2.sort()
+ b3 = bo3["index","morton_index"]
+ b3.sort()
+ assert_array_equal(b1, overlap)
+ assert_array_equal(b2, diff)
+ assert_array_equal(b3, both)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ assert_array_equal(b3, b4)
+ assert_array_equal(b1, b5)
+ bo6 = ell1 ^ ell2
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(i1, i2))
+
+def test_boolean_mix_periodicity():
+ r"""Test that a hybrid boolean region behaves as we expect.
+
+ This also tests nested logic and that periodicity works.
+ """
+ ds = fake_amr_ds()
+ re = ds.region([0.5]*3, [0.0]*3, [1]*3) # whole thing
+ sp = ds.sphere([0.95]*3, 0.3) # wraps around
+ cyl = ds.disk([0.05]*3, [1,1,1], 0.1, 0.4) # wraps around
+ # Get original indices
+ rei = re["index","morton_index"]
+ spi = sp["index","morton_index"]
+ cyli = cyl["index","morton_index"]
+ # Make some booleans
+ # whole box minux spherical bites at corners
+ bo1 = re - sp
+ # sphere plus cylinder
+ bo2 = sp | cyl
+ # a jumble, the region minus the sp+cyl
+ bo3 = re - (sp | cyl)
+ # Now make sure the indices also behave as we expect.
+ bo4 = ds.union([re, sp, cyl])
+ bo5 = ds.intersection([re, sp, cyl])
+ expect = np.setdiff1d(rei, spi)
+ ii = bo1["index","morton_index"]
+ ii.sort()
+ assert_array_equal(expect, ii)
+ #
+ expect = np.union1d(spi, cyli)
+ ii = bo2["index","morton_index"]
+ ii.sort()
+ assert_array_equal(expect, ii)
+ #
+ expect = np.union1d(spi, cyli)
+ expect = np.setdiff1d(rei, expect)
+ ii = bo3["index","morton_index"]
+ ii.sort()
+ assert_array_equal(expect, ii)
+ b4 = bo4["index","morton_index"]
+ b4.sort()
+ b5 = bo5["index","morton_index"]
+ b5.sort()
+ ii = np.union1d(np.union1d(rei, cyli), spi)
+ ii.sort()
+ assert_array_equal(ii, b4)
+ ii = np.intersect1d(np.intersect1d(rei, cyli), spi)
+ ii.sort()
+ assert_array_equal(ii, b5)
+
+ bo6 = (re ^ sp) ^ cyl
+ b6 = bo6["index", "morton_index"]
+ b6.sort()
+ assert_array_equal(b6, np.setxor1d(np.setxor1d(rei, spi), cyli))
+
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/angular_momentum.py
--- a/yt/fields/angular_momentum.py
+++ b/yt/fields/angular_momentum.py
@@ -59,15 +59,15 @@
rv = data.ds.arr(rv, input_units = data["index", "x"].units)
return xv * rv[...,1] - yv * rv[...,0]
- registry.add_field((ftype, "specific_angular_momentum_x"),
+ registry.add_field((ftype, "specific_angular_momentum_x"), sampling_type="cell",
function=_specific_angular_momentum_x,
units=unit_system["specific_angular_momentum"],
validators=[ValidateParameter("center")])
- registry.add_field((ftype, "specific_angular_momentum_y"),
+ registry.add_field((ftype, "specific_angular_momentum_y"), sampling_type="cell",
function=_specific_angular_momentum_y,
units=unit_system["specific_angular_momentum"],
validators=[ValidateParameter("center")])
- registry.add_field((ftype, "specific_angular_momentum_z"),
+ registry.add_field((ftype, "specific_angular_momentum_z"), sampling_type="cell",
function=_specific_angular_momentum_z,
units=unit_system["specific_angular_momentum"],
validators=[ValidateParameter("center")])
@@ -78,7 +78,7 @@
def _angular_momentum_x(field, data):
return data[ftype, "cell_mass"] \
* data[ftype, "specific_angular_momentum_x"]
- registry.add_field((ftype, "angular_momentum_x"),
+ registry.add_field((ftype, "angular_momentum_x"), sampling_type="cell",
function=_angular_momentum_x,
units=unit_system["angular_momentum"],
validators=[ValidateParameter('center')])
@@ -86,7 +86,7 @@
def _angular_momentum_y(field, data):
return data[ftype, "cell_mass"] \
* data[ftype, "specific_angular_momentum_y"]
- registry.add_field((ftype, "angular_momentum_y"),
+ registry.add_field((ftype, "angular_momentum_y"), sampling_type="cell",
function=_angular_momentum_y,
units=unit_system["angular_momentum"],
validators=[ValidateParameter('center')])
@@ -94,7 +94,7 @@
def _angular_momentum_z(field, data):
return data[ftype, "cell_mass"] \
* data[ftype, "specific_angular_momentum_z"]
- registry.add_field((ftype, "angular_momentum_z"),
+ registry.add_field((ftype, "angular_momentum_z"), sampling_type="cell",
function=_angular_momentum_z,
units=unit_system["angular_momentum"],
validators=[ValidateParameter('center')])
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/astro_fields.py
--- a/yt/fields/astro_fields.py
+++ b/yt/fields/astro_fields.py
@@ -52,7 +52,7 @@
"""
return np.sqrt(3.0 * np.pi / (16.0 * G * data[ftype, "density"]))
- registry.add_field((ftype, "dynamical_time"),
+ registry.add_field((ftype, "dynamical_time"), sampling_type="cell",
function=_dynamical_time,
units=unit_system["time"])
@@ -65,7 +65,7 @@
(data[ftype, "density"]**(-0.5)))
return u
- registry.add_field((ftype, "jeans_mass"),
+ registry.add_field((ftype, "jeans_mass"), sampling_type="cell",
function=_jeans_mass,
units=unit_system["mass"])
@@ -88,7 +88,7 @@
- 1.6667 * logT0**1 - 0.2193 * logT0)),
"") # add correct units here
- registry.add_field((ftype, "chandra_emissivity"),
+ registry.add_field((ftype, "chandra_emissivity"), sampling_type="cell",
function=_chandra_emissivity,
units="") # add correct units here
@@ -102,7 +102,7 @@
nenh *= 0.5*(1.+X_H)*X_H*data["cell_volume"]
return nenh
- registry.add_field((ftype, "emission_measure"),
+ registry.add_field((ftype, "emission_measure"), sampling_type="cell",
function=_emission_measure,
units=unit_system["number_density"])
@@ -112,7 +112,7 @@
* data[ftype, "temperature"].to_ndarray()**0.5,
"") # add correct units here
- registry.add_field((ftype, "xray_emissivity"),
+ registry.add_field((ftype, "xray_emissivity"), sampling_type="cell",
function=_xray_emissivity,
units="") # add correct units here
@@ -121,7 +121,7 @@
# Only useful as a weight_field for temperature, metallicity, velocity
return data["density"]*data["density"]*data["kT"]**-0.25/mh/mh
- registry.add_field((ftype,"mazzotta_weighting"),
+ registry.add_field((ftype,"mazzotta_weighting"), sampling_type="cell",
function=_mazzotta_weighting,
units="keV**-0.25*cm**-6")
@@ -135,7 +135,7 @@
# See issue #1225
return -scale * vel * data[ftype, "density"]
- registry.add_field((ftype, "sz_kinetic"),
+ registry.add_field((ftype, "sz_kinetic"), sampling_type="cell",
function=_sz_kinetic,
units=unit_system["length"]**-1,
validators=[
@@ -145,6 +145,6 @@
scale = 0.88 / mh * kboltz / (me * clight*clight) * sigma_thompson
return scale * data[ftype, "density"] * data[ftype, "temperature"]
- registry.add_field((ftype, "szy"),
+ registry.add_field((ftype, "szy"), sampling_type="cell",
function=_szy,
units=unit_system["length"]**-1)
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/cosmology_fields.py
--- a/yt/fields/cosmology_fields.py
+++ b/yt/fields/cosmology_fields.py
@@ -45,14 +45,14 @@
return data[ftype, "density"] + \
data[ftype, "dark_matter_density"]
- registry.add_field((ftype, "matter_density"),
+ registry.add_field((ftype, "matter_density"), sampling_type="cell",
function=_matter_density,
units=unit_system["density"])
def _matter_mass(field, data):
return data[ftype, "matter_density"] * data["index", "cell_volume"]
- registry.add_field((ftype, "matter_mass"),
+ registry.add_field((ftype, "matter_mass"), sampling_type="cell",
function=_matter_mass,
units=unit_system["mass"])
@@ -65,7 +65,7 @@
return data[ftype, "matter_density"] / \
co.critical_density(data.ds.current_redshift)
- registry.add_field((ftype, "overdensity"),
+ registry.add_field((ftype, "overdensity"), sampling_type="cell",
function=_overdensity,
units="")
@@ -83,7 +83,7 @@
return data[ftype, "density"] / omega_baryon / co.critical_density(0.0) / \
(1.0 + data.ds.current_redshift)**3
- registry.add_field((ftype, "baryon_overdensity"),
+ registry.add_field((ftype, "baryon_overdensity"), sampling_type="cell",
function=_baryon_overdensity,
units="",
validators=[ValidateParameter("omega_baryon")])
@@ -100,7 +100,7 @@
co.critical_density(0.0) / \
(1.0 + data.ds.current_redshift)**3
- registry.add_field((ftype, "matter_overdensity"),
+ registry.add_field((ftype, "matter_overdensity"), sampling_type="cell",
function=_matter_overdensity,
units="")
@@ -109,7 +109,7 @@
virial_radius = data.get_field_parameter("virial_radius")
return data["radius"] / virial_radius
- registry.add_field(("index", "virial_radius_fraction"),
+ registry.add_field(("index", "virial_radius_fraction"), sampling_type="cell",
function=_virial_radius_fraction,
validators=[ValidateParameter("virial_radius")],
units="")
@@ -137,7 +137,7 @@
return (1.5 * (co.hubble_constant / speed_of_light_cgs)**2 * (dl * dls / ds) * \
data[ftype, "matter_overdensity"]).in_units("1/cm")
- registry.add_field((ftype, "weak_lensing_convergence"),
+ registry.add_field((ftype, "weak_lensing_convergence"), sampling_type="cell",
function=_weak_lensing_convergence,
units=unit_system["length"]**-1,
validators=[ValidateParameter("observer_redshift"),
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/derived_field.py
--- a/yt/fields/derived_field.py
+++ b/yt/fields/derived_field.py
@@ -91,8 +91,8 @@
The dimensions of the field, only needed if units="auto" and only used
for error checking.
"""
- def __init__(self, name, function, units=None,
- take_log=True, validators=None, sampling_type = "mesh",
+ def __init__(self, name, sampling_type, function, units=None,
+ take_log=True, validators=None,
particle_type=None, vector_field=False, display_field=True,
not_in_all=False, display_name=None, output_units=None,
dimensions=None, ds=None):
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -97,9 +97,9 @@
output_units = units
if (ptype, f) not in self.field_list:
continue
- self.add_output_field((ptype, f),
- units = units, sampling_type = "particle",
- display_name = dn, output_units = output_units)
+ self.add_output_field((ptype, f), sampling_type="particle",
+ units = units, display_name = dn,
+ output_units = output_units)
for alias in aliases:
self.alias((ptype, alias), (ptype, f), units = output_units)
@@ -134,9 +134,8 @@
raise RuntimeError
if field[0] not in self.ds.particle_types:
continue
- self.add_output_field(field,
- units = self.ds.field_units.get(field, ""),
- sampling_type = "particle")
+ self.add_output_field(field, sampling_type="particle",
+ units = self.ds.field_units.get(field, ""))
self.setup_smoothed_fields(ptype,
num_neighbors=num_neighbors,
ftype=ftype)
@@ -196,12 +195,12 @@
units = ""
elif units == 1.0:
units = ""
- self.add_output_field(field, units = units,
+ self.add_output_field(field, sampling_type="cell", units = units,
display_name = display_name)
for alias in aliases:
self.alias(("gas", alias), field)
- def add_field(self, name, function=None, **kwargs):
+ def add_field(self, name, sampling_type, function=None, **kwargs):
"""
Add a new field, along with supplemental metadata, to the list of
available fields. This respects a number of arguments, all of which
@@ -250,12 +249,12 @@
kwargs.setdefault('ds', self.ds)
if function is None:
def create_function(f):
- self[name] = DerivedField(name, f, **kwargs)
+ self[name] = DerivedField(name, sampling_type, f, **kwargs)
return f
return create_function
if isinstance(name, tuple):
- self[name] = DerivedField(name, function, **kwargs)
+ self[name] = DerivedField(name, sampling_type, function, **kwargs)
return
if kwargs.get("particle_type", False):
@@ -265,10 +264,11 @@
if (ftype, name) not in self:
tuple_name = (ftype, name)
- self[tuple_name] = DerivedField(tuple_name, function, **kwargs)
+ self[tuple_name] = DerivedField(tuple_name, sampling_type, function,
+ **kwargs)
self.alias(name, tuple_name)
else:
- self[name] = DerivedField(name, function, **kwargs)
+ self[name] = DerivedField(name, sampling_type, function, **kwargs)
def load_all_plugins(self, ftype="gas"):
loaded = []
@@ -296,9 +296,9 @@
self.ds.derived_field_list = list(sorted(dfl, key=tupleize))
return loaded, unavailable
- def add_output_field(self, name, **kwargs):
+ def add_output_field(self, name, sampling_type, **kwargs):
kwargs.setdefault('ds', self.ds)
- self[name] = DerivedField(name, NullFunc, **kwargs)
+ self[name] = DerivedField(name, sampling_type, NullFunc, **kwargs)
def alias(self, alias_name, original_name, units = None):
if original_name not in self: return
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/fluid_fields.py
--- a/yt/fields/fluid_fields.py
+++ b/yt/fields/fluid_fields.py
@@ -59,14 +59,14 @@
def _cell_mass(field, data):
return data[ftype, "density"] * data[ftype, "cell_volume"]
- registry.add_field((ftype, "cell_mass"),
+ registry.add_field((ftype, "cell_mass"), sampling_type="cell",
function=_cell_mass,
units=unit_system["mass"])
def _sound_speed(field, data):
tr = data.ds.gamma * data[ftype, "pressure"] / data[ftype, "density"]
return np.sqrt(tr)
- registry.add_field((ftype, "sound_speed"),
+ registry.add_field((ftype, "sound_speed"), sampling_type="cell",
function=_sound_speed,
units=unit_system["velocity"])
@@ -74,7 +74,7 @@
""" Radial component of M{|v|/c_sound} """
tr = data[ftype, "radial_velocity"] / data[ftype, "sound_speed"]
return np.abs(tr)
- registry.add_field((ftype, "radial_mach_number"),
+ registry.add_field((ftype, "radial_mach_number"), sampling_type="cell",
function=_radial_mach_number,
units = "")
@@ -82,14 +82,14 @@
return 0.5*data[ftype, "density"] * ( data[ftype, "velocity_x"]**2.0
+ data[ftype, "velocity_y"]**2.0
+ data[ftype, "velocity_z"]**2.0 )
- registry.add_field((ftype, "kinetic_energy"),
+ registry.add_field((ftype, "kinetic_energy"), sampling_type="cell",
function = _kin_energy,
units = unit_system["pressure"])
def _mach_number(field, data):
""" M{|v|/c_sound} """
return data[ftype, "velocity_magnitude"] / data[ftype, "sound_speed"]
- registry.add_field((ftype, "mach_number"),
+ registry.add_field((ftype, "mach_number"), sampling_type="cell",
function=_mach_number,
units = "")
@@ -103,7 +103,7 @@
tr = np.minimum(np.minimum(t1, t2), t3)
return tr
- registry.add_field((ftype, "courant_time_step"),
+ registry.add_field((ftype, "courant_time_step"), sampling_type="cell",
function=_courant_time_step,
units=unit_system["time"])
@@ -113,13 +113,13 @@
* (data[ftype, "density"] * data[ftype, "thermal_energy"])
return tr
- registry.add_field((ftype, "pressure"),
+ registry.add_field((ftype, "pressure"), sampling_type="cell",
function=_pressure,
units=unit_system["pressure"])
def _kT(field, data):
return (kboltz*data[ftype, "temperature"]).in_units("keV")
- registry.add_field((ftype, "kT"),
+ registry.add_field((ftype, "kT"), sampling_type="cell",
function=_kT,
units="keV",
display_name="Temperature")
@@ -132,7 +132,7 @@
gammam1 = 2./3.
tr = data[ftype,"kT"] / ((data[ftype, "density"]/mw)**gammam1)
return data.apply_units(tr, field.units)
- registry.add_field((ftype, "entropy"),
+ registry.add_field((ftype, "entropy"), sampling_type="cell",
units="keV*cm**2",
function=_entropy)
@@ -140,13 +140,13 @@
tr = data[ftype, "metal_density"] / data[ftype, "density"]
tr /= metallicity_sun
return data.apply_units(tr, "Zsun")
- registry.add_field((ftype, "metallicity"),
+ registry.add_field((ftype, "metallicity"), sampling_type="cell",
function=_metallicity,
units="Zsun")
def _metal_mass(field, data):
return data[ftype, "metal_density"] * data[ftype, "cell_volume"]
- registry.add_field((ftype, "metal_mass"),
+ registry.add_field((ftype, "metal_mass"), sampling_type="cell",
function=_metal_mass,
units=unit_system["mass"])
@@ -156,13 +156,13 @@
for species in data.ds.field_info.species_names:
field_data += data["gas", "%s_number_density" % species]
return field_data
- registry.add_field((ftype, "number_density"),
+ registry.add_field((ftype, "number_density"), sampling_type="cell",
function = _number_density,
units=unit_system["number_density"])
def _mean_molecular_weight(field, data):
return (data[ftype, "density"] / (mh * data[ftype, "number_density"]))
- registry.add_field((ftype, "mean_molecular_weight"),
+ registry.add_field((ftype, "mean_molecular_weight"), sampling_type="cell",
function=_mean_molecular_weight,
units="")
@@ -207,7 +207,7 @@
for axi, ax in enumerate('xyz'):
f = grad_func(axi, ax)
- registry.add_field((ftype, "%s_gradient_%s" % (fname, ax)),
+ registry.add_field((ftype, "%s_gradient_%s" % (fname, ax)), sampling_type="cell",
function = f,
validators = [ValidateSpatial(1, [grad_field])],
units = grad_units)
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/fluid_vector_fields.py
--- a/yt/fields/fluid_vector_fields.py
+++ b/yt/fields/fluid_vector_fields.py
@@ -67,7 +67,7 @@
bv_validators = [ValidateSpatial(1, [(ftype, "density"), (ftype, "pressure")])]
for ax in 'xyz':
n = "baroclinic_vorticity_%s" % ax
- registry.add_field((ftype, n), function=eval("_%s" % n),
+ registry.add_field((ftype, n), sampling_type="cell", function=eval("_%s" % n),
validators=bv_validators,
units=unit_system["frequency"]**2)
@@ -119,7 +119,7 @@
(ftype, "velocity_z")])]
for ax in 'xyz':
n = "vorticity_%s" % ax
- registry.add_field((ftype, n),
+ registry.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units=unit_system["frequency"],
validators=vort_validators)
@@ -138,7 +138,7 @@
return data[ftype, "velocity_divergence"] * data[ftype, "vorticity_z"]
for ax in 'xyz':
n = "vorticity_stretching_%s" % ax
- registry.add_field((ftype, n),
+ registry.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units = unit_system["frequency"]**2,
validators=vort_validators)
@@ -159,7 +159,7 @@
data[ftype, "baroclinic_vorticity_z"]
for ax in 'xyz':
n = "vorticity_growth_%s" % ax
- registry.add_field((ftype, n),
+ registry.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units=unit_system["frequency"]**2,
validators=vort_validators)
@@ -175,7 +175,7 @@
result = np.sign(dot) * result
return result
- registry.add_field((ftype, "vorticity_growth_magnitude"),
+ registry.add_field((ftype, "vorticity_growth_magnitude"), sampling_type="cell",
function=_vorticity_growth_magnitude,
units=unit_system["frequency"]**2,
validators=vort_validators,
@@ -186,7 +186,7 @@
data[ftype, "vorticity_growth_y"]**2 +
data[ftype, "vorticity_growth_z"]**2)
- registry.add_field((ftype, "vorticity_growth_magnitude_absolute"),
+ registry.add_field((ftype, "vorticity_growth_magnitude_absolute"), sampling_type="cell",
function=_vorticity_growth_magnitude_absolute,
units=unit_system["frequency"]**2,
validators=vort_validators)
@@ -197,7 +197,7 @@
domegaz_dt = data[ftype, "vorticity_z"] / data[ftype, "vorticity_growth_z"]
return np.sqrt(domegax_dt**2 + domegay_dt**2 + domegaz_dt**2)
- registry.add_field((ftype, "vorticity_growth_timescale"),
+ registry.add_field((ftype, "vorticity_growth_timescale"), sampling_type="cell",
function=_vorticity_growth_timescale,
units=unit_system["time"],
validators=vort_validators)
@@ -232,7 +232,7 @@
(ftype, "radiation_acceleration_z")])]
for ax in 'xyz':
n = "vorticity_radiation_pressure_%s" % ax
- registry.add_field((ftype, n),
+ registry.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units=unit_system["frequency"]**2,
validators=vrp_validators)
@@ -257,7 +257,7 @@
for ax in 'xyz':
n = "vorticity_radiation_pressure_growth_%s" % ax
- registry.add_field((ftype, n),
+ registry.add_field((ftype, n), sampling_type="cell",
function=eval("_%s" % n),
units=unit_system["frequency"]**2,
validators=vrp_validators)
@@ -273,7 +273,7 @@
result = np.sign(dot) * result
return result
- registry.add_field((ftype, "vorticity_radiation_pressure_growth_magnitude"),
+ registry.add_field((ftype, "vorticity_radiation_pressure_growth_magnitude"), sampling_type="cell",
function=_vorticity_radiation_pressure_growth_magnitude,
units=unit_system["frequency"]**2,
validators=vrp_validators,
@@ -284,7 +284,7 @@
data[ftype, "vorticity_radiation_pressure_growth_y"]**2 +
data[ftype, "vorticity_radiation_pressure_growth_z"]**2)
- registry.add_field((ftype, "vorticity_radiation_pressure_growth_magnitude_absolute"),
+ registry.add_field((ftype, "vorticity_radiation_pressure_growth_magnitude_absolute"), sampling_type="cell",
function=_vorticity_radiation_pressure_growth_magnitude_absolute,
units="s**(-2)",
validators=vrp_validators)
@@ -298,7 +298,7 @@
data[ftype, "vorticity_radiation_pressure_growth_z"]
return np.sqrt(domegax_dt**2 + domegay_dt**2 + domegaz_dt**2)
- registry.add_field((ftype, "vorticity_radiation_pressure_growth_timescale"),
+ registry.add_field((ftype, "vorticity_radiation_pressure_growth_timescale"), sampling_type="cell",
function=_vorticity_radiation_pressure_growth_timescale,
units=unit_system["time"],
validators=vrp_validators)
@@ -344,7 +344,7 @@
new_field[sl_center, sl_center, sl_center] = f
return new_field
- registry.add_field((ftype, "shear"),
+ registry.add_field((ftype, "shear"), sampling_type="cell",
function=_shear,
validators=[ValidateSpatial(1,
[(ftype, "velocity_x"),
@@ -361,7 +361,7 @@
return data[ftype, "shear"] / data[ftype, "sound_speed"]
- registry.add_field((ftype, "shear_criterion"),
+ registry.add_field((ftype, "shear_criterion"), sampling_type="cell",
function=_shear_criterion,
units=unit_system["length"]**-1,
validators=[ValidateSpatial(1,
@@ -417,7 +417,7 @@
new_field[sl_center, sl_center, sl_center] = f
return new_field
- registry.add_field((ftype, "shear_mach"),
+ registry.add_field((ftype, "shear_mach"), sampling_type="cell",
function=_shear_mach,
units="",
validators=[ValidateSpatial(1,
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/geometric_fields.py
--- a/yt/fields/geometric_fields.py
+++ b/yt/fields/geometric_fields.py
@@ -46,7 +46,7 @@
"""
return get_radius(data, "")
- registry.add_field(("index", "radius"),
+ registry.add_field(("index", "radius"), sampling_type="cell",
function=_radius,
validators=[ValidateParameter("center")],
units=unit_system["length"])
@@ -59,7 +59,7 @@
return data._reshape_vals(arr)
return arr
- registry.add_field(("index", "grid_level"),
+ registry.add_field(("index", "grid_level"), sampling_type="cell",
function=_grid_level,
units="",
validators=[ValidateSpatial(0)])
@@ -76,7 +76,7 @@
return data._reshape_vals(arr)
return arr
- registry.add_field(("index", "grid_indices"),
+ registry.add_field(("index", "grid_indices"), sampling_type="cell",
function=_grid_indices,
units="",
validators=[ValidateSpatial(0)],
@@ -87,7 +87,7 @@
return np.ones(data["index", "ones"].shape,
dtype="float64")/data["index", "dx"]
- registry.add_field(("index", "ones_over_dx"),
+ registry.add_field(("index", "ones_over_dx"), sampling_type="cell",
function=_ones_over_dx,
units=unit_system["length"]**-1,
display_field=False)
@@ -97,7 +97,7 @@
arr = np.zeros(data["index", "ones"].shape, dtype='float64')
return data.apply_units(arr, field.units)
- registry.add_field(("index", "zeros"),
+ registry.add_field(("index", "zeros"), sampling_type="cell",
function=_zeros,
units="",
display_field=False)
@@ -109,7 +109,7 @@
return data._reshape_vals(arr)
return data.apply_units(arr, field.units)
- registry.add_field(("index", "ones"),
+ registry.add_field(("index", "ones"), sampling_type="cell",
function=_ones,
units="",
display_field=False)
@@ -132,7 +132,7 @@
data["index", "z"].ravel(), LE, RE)
morton.shape = data["index", "x"].shape
return morton.view("f8")
- registry.add_field(("index", "morton_index"), function=_morton_index,
+ registry.add_field(("index", "morton_index"), sampling_type="cell", function=_morton_index,
units = "")
def _spherical_radius(field, data):
@@ -144,7 +144,7 @@
coords = get_periodic_rvec(data)
return data.ds.arr(get_sph_r(coords), "code_length").in_base(unit_system.name)
- registry.add_field(("index", "spherical_radius"),
+ registry.add_field(("index", "spherical_radius"), sampling_type="cell",
function=_spherical_radius,
validators=[ValidateParameter("center")],
units=unit_system["length"])
@@ -153,7 +153,7 @@
"""This field is deprecated and will be removed in a future release"""
return data['index', 'spherical_radius']
- registry.add_field(("index", "spherical_r"),
+ registry.add_field(("index", "spherical_r"), sampling_type="cell",
function=_spherical_r,
validators=[ValidateParameter("center")],
units=unit_system["length"])
@@ -171,7 +171,7 @@
coords = get_periodic_rvec(data)
return get_sph_theta(coords, normal)
- registry.add_field(("index", "spherical_theta"),
+ registry.add_field(("index", "spherical_theta"), sampling_type="cell",
function=_spherical_theta,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
@@ -190,7 +190,7 @@
coords = get_periodic_rvec(data)
return get_sph_phi(coords, normal)
- registry.add_field(("index", "spherical_phi"),
+ registry.add_field(("index", "spherical_phi"), sampling_type="cell",
function=_spherical_phi,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
@@ -206,7 +206,7 @@
coords = get_periodic_rvec(data)
return data.ds.arr(get_cyl_r(coords, normal), "code_length").in_base(unit_system.name)
- registry.add_field(("index", "cylindrical_radius"),
+ registry.add_field(("index", "cylindrical_radius"), sampling_type="cell",
function=_cylindrical_radius,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
@@ -216,7 +216,7 @@
"""This field is deprecated and will be removed in a future release"""
return data['index', 'cylindrical_radius']
- registry.add_field(("index", "cylindrical_r"),
+ registry.add_field(("index", "cylindrical_r"), sampling_type="cell",
function=_cylindrical_r,
validators=[ValidateParameter("center")],
units=unit_system["length"])
@@ -231,7 +231,7 @@
coords = get_periodic_rvec(data)
return data.ds.arr(get_cyl_z(coords, normal), "code_length").in_base(unit_system.name)
- registry.add_field(("index", "cylindrical_z"),
+ registry.add_field(("index", "cylindrical_z"), sampling_type="cell",
function=_cylindrical_z,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
@@ -250,7 +250,7 @@
coords = get_periodic_rvec(data)
return get_cyl_theta(coords, normal)
- registry.add_field(("index", "cylindrical_theta"),
+ registry.add_field(("index", "cylindrical_theta"), sampling_type="cell",
function=_cylindrical_theta,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
@@ -260,7 +260,7 @@
"""This field is dprecated and will be removed in a future release"""
return data["index", "spherical_theta"]
- registry.add_field(("index", "disk_angle"),
+ registry.add_field(("index", "disk_angle"), sampling_type="cell",
function=_disk_angle,
take_log=False,
display_field=False,
@@ -272,7 +272,7 @@
"""This field is deprecated and will be removed in a future release"""
return data["index", "cylindrical_z"]
- registry.add_field(("index", "height"),
+ registry.add_field(("index", "height"), sampling_type="cell",
function=_height,
validators=[ValidateParameter("center"),
ValidateParameter("normal")],
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/local_fields.py
--- a/yt/fields/local_fields.py
+++ b/yt/fields/local_fields.py
@@ -12,6 +12,7 @@
#
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+import warnings
from yt.utilities.logger import \
ytLogger as mylog
@@ -23,7 +24,7 @@
FieldInfoContainer
class LocalFieldInfoContainer(FieldInfoContainer):
- def add_field(self, name, function=None, **kwargs):
+ def add_field(self, name, function=None, sampling_type=None, **kwargs):
if not isinstance(name, tuple):
if kwargs.setdefault('particle_type', False):
name = ('all', name)
@@ -34,8 +35,19 @@
if not override and name in self:
mylog.warning("Field %s already exists. To override use " +
"force_override=True.", name)
+ if kwargs.setdefault('particle_type', False):
+ if sampling_type is not None and sampling_type != "particle":
+ raise RuntimeError("Clashing definition of 'sampling_type' and "
+ "'particle_type'. Note that 'particle_type' is "
+ "deprecated. Please just use 'sampling_type'.")
+ else:
+ sampling_type = "particle"
+ if sampling_type is None:
+ warnings.warn("Because 'sampling_type' not specified, yt will "
+ "assume a cell 'sampling_type'")
+ sampling_type = "cell"
return super(LocalFieldInfoContainer,
- self).add_field(name, function, **kwargs)
+ self).add_field(name, sampling_type, function, **kwargs)
# Empty FieldInfoContainer
local_fields = LocalFieldInfoContainer(None, [], None)
diff -r 4f3f870ced4e4d5785c9a7b2562b9c22f563f91e -r cec7d20114951c59abdaa02f7527a5db19b93434 yt/fields/magnetic_field.py
--- a/yt/fields/magnetic_field.py
+++ b/yt/fields/magnetic_field.py
@@ -46,26 +46,26 @@
data[ftype,"magnetic_field_y"]**2 +
data[ftype,"magnetic_field_z"]**2)
return np.sqrt(B2)
- registry.add_field((ftype,"magnetic_field_strength"),
+ registry.add_field((ftype,"magnetic_field_strength"), sampling_type="cell",
function=_magnetic_field_strength,
units=u)
def _magnetic_energy(field, data):
B = data[ftype,"magnetic_field_strength"]
return 0.5*B*B/mag_factors[B.units.dimensions]
- registry.add_field((ftype, "magnetic_energy"),
+ registry.add_field((ftype, "magnetic_energy"), sampling_type="cell",
function=_magnetic_energy,
units=unit_system["pressure"])
def _plasma_beta(field,data):
return data[ftype,'pressure']/data[ftype,'magnetic_energy']
- registry.add_field((ftype, "plasma_beta"),
+ registry.add_field((ftype, "plasma_beta"), sampling_type="cell",
function=_plasma_beta,
units="")
def _magnetic_pressure(field,data):
return data[ftype,'magnetic_energy']
- registry.add_field((ftype, "magnetic_pressure"),
+ registry.add_field((ftype, "magnetic_pressure"), sampling_type="cell",
function=_magnetic_pressure,
units=unit_system["pressure"])
@@ -83,7 +83,7 @@
return get_sph_theta_component(Bfields, theta, phi, normal)
- registry.add_field((ftype, "magnetic_field_poloidal"),
+ registry.add_field((ftype, "magnetic_field_poloidal"), sampling_type="cell",
function=_magnetic_field_poloidal,
units=u, validators=[ValidateParameter("normal")])
@@ -99,19 +99,19 @@
phi = data["index", 'spherical_phi']
return get_sph_phi_component(Bfields, phi, normal)
- registry.add_field((ftype, "magnetic_field_toroidal"),
+ registry.add_field((ftype, "magnetic_field_toroidal"), sampling_type="cell",
function=_magnetic_field_toroidal,
units=u, validators=[ValidateParameter("normal")])
def _alfven_speed(field,data):
B = data[ftype,'magnetic_field_strength']
return B/np.sqrt(mag_factors[B.units.dimensions]*data[ftype,'density'])
- registry.add_field((ftype, "alfven_speed"), function=_alfven_speed,
+ registry.add_field((ftype, "alfven_speed"), sampling_type="cell", function=_alfven_speed,
units=unit_system["velocity"])
def _mach_alfven(field,data):
return data[ftype,'velocity_magnitude']/data[ftype,'alfven_speed']
- registry.add_field((ftype, "mach_alfven"), function=_mach_alfven,
+ registry.add_field((ftype, "mach_alfven"), sampling_type="cell", function=_mach_alfven,
units="dimensionless")
def setup_magnetic_field_aliases(registry, ds_ftype, ds_fields, ftype="gas"):
@@ -161,6 +161,6 @@
return convert(data[fd])
return _mag_field
for ax, fd in zip("xyz", ds_fields):
- registry.add_field((ftype,"magnetic_field_%s" % ax),
+ registry.add_field((ftype,"magnetic_field_%s" % ax), sampling_type="cell",
function=mag_field(fd),
- units=unit_system[to_units.dimensions])
\ No newline at end of file
+ units=unit_system[to_units.dimensions])
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/be8a962c34c6/
Changeset: be8a962c34c6
Branch: yt
User: al007
Date: 2016-09-30 21:32:08+00:00
Summary: Merge north_vector.
Affected #: 2 files
diff -r cec7d20114951c59abdaa02f7527a5db19b93434 -r be8a962c34c60755fabbe1d6add6ec54e195d4f1 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -179,7 +179,7 @@
'Cannot set the camera focus and position to the same value')
self._position = position
self.switch_orientation(normal_vector=self.focus - self._position,
- north_vector=None)
+ north_vector=self.north_vector)
def fdel(self):
del self._position
@@ -386,9 +386,9 @@
calculated automatically.
"""
+ if north_vector is not None:
+ self.north_vector = north_vector
self.position = position
- self.switch_orientation(normal_vector=self.focus - self.position,
- north_vector=north_vector)
def get_position(self):
"""Return the current camera position"""
diff -r cec7d20114951c59abdaa02f7527a5db19b93434 -r be8a962c34c60755fabbe1d6add6ec54e195d4f1 yt/visualization/volume_rendering/tests/test_camera_attributes.py
--- a/yt/visualization/volume_rendering/tests/test_camera_attributes.py
+++ b/yt/visualization/volume_rendering/tests/test_camera_attributes.py
@@ -110,3 +110,12 @@
for lens_type in valid_lens_types:
cam.set_lens(lens_type)
+
+ # See issue #1287
+ cam.focus = [0, 0, 0]
+ cam_pos = [1, 0, 0]
+ north_vector = [0, 1, 0]
+ cam.set_position(cam_pos, north_vector)
+ cam_pos = [0, 1, 0]
+ north_vector = [0, 0, 1]
+ cam.set_position(cam_pos, north_vector)
https://bitbucket.org/yt_analysis/yt/commits/f02e8a513aa7/
Changeset: f02e8a513aa7
Branch: yt
User: al007
Date: 2016-09-30 22:11:14+00:00
Summary: Remove commented out tet10 lines in render_source.py for readability.
Affected #: 1 file
diff -r be8a962c34c60755fabbe1d6add6ec54e195d4f1 -r f02e8a513aa7b68f5f195001a89a2f9c6669d0f5 yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -620,12 +620,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- # elif indices.shape[1] == 10:
- # # tetrahedral
- # mylog.warning("10-node tetrahedral elements not yet supported, " +
- # "dropping to 1st order.")
- # field_data = field_data[:, 0:4]
- # indices = indices[:, 0:4]
self.mesh = mesh_construction.LinearElementMesh(self.volume,
vertices,
@@ -659,12 +653,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- # elif indices.shape[1] == 10:
- # # tetrahedral
- # mylog.warning("10-node tetrahedral elements not yet supported, " +
- # "dropping to 1st order.")
- # field_data = field_data[:, 0:4]
- # indices = indices[:, 0:4]
self.volume = BVH(vertices, indices, field_data)
https://bitbucket.org/yt_analysis/yt/commits/897ffcfa02ef/
Changeset: 897ffcfa02ef
Branch: yt
User: al007
Date: 2016-10-03 22:31:10+00:00
Summary: Fix u,v bounds so ray hits correctly ID'd.
Affected #: 3 files
diff -r f02e8a513aa7b68f5f195001a89a2f9c6669d0f5 -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -834,11 +834,10 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tris, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(2):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] - 1.0 > self.inclusion_tol):
+ return 0
return 1
cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
@@ -893,11 +892,12 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tets, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(3):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[2] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] + mapped_coord[2] - 1.0 > \
+ self.inclusion_tol):
+ return 0
return 1
@cython.boundscheck(False)
diff -r f02e8a513aa7b68f5f195001a89a2f9c6669d0f5 -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe yt/utilities/lib/mesh_triangulation.pyx
--- a/yt/utilities/lib/mesh_triangulation.pyx
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -180,7 +180,7 @@
if (self.VPE == 8 or self.VPE == 20 or self.VPE == 27):
self.TPE = HEX_NT
self.tri_array = triangulate_hex
- elif self.VPE == 4:
+ elif (self.VPE == 4 or self.VPE == 10):
self.TPE = TETRA_NT
self.tri_array = triangulate_tetra
elif self.VPE == 6:
diff -r f02e8a513aa7b68f5f195001a89a2f9c6669d0f5 -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -445,7 +445,7 @@
if (hd.t < ray.t_near or hd.t > ray.t_far):
return False
- if ((0 <= hd.u <= 1.0) and (0 <= hd.v <= 1.0) and hd.converged):
+ if (0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1.0 and hd.converged):
# we have a hit, so update ray information
ray.t_far = hd.t
ray.elem_id = tet_patch.elem_id
https://bitbucket.org/yt_analysis/yt/commits/f42ee7dc1a88/
Changeset: f42ee7dc1a88
Branch: yt
User: al007
Date: 2016-10-04 23:00:03+00:00
Summary: Adding embree support.
Affected #: 6 files
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_construction.pxd
--- a/yt/utilities/lib/mesh_construction.pxd
+++ b/yt/utilities/lib/mesh_construction.pxd
@@ -18,3 +18,9 @@
unsigned int geomID
np.float64_t* vertices
np.float64_t* field_data
+
+ctypedef struct Tet_Patch:
+ float[6][3] v
+ unsigned int geomID
+ np.float64_t* vertices
+ np.float64_t* field_data
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -37,7 +37,9 @@
sample_wedge
from mesh_intersection cimport \
patchIntersectFunc, \
- patchBoundsFunc
+ patchBoundsFunc, \
+ tet_patchIntersectFunc, \
+ tet_patchBoundsFunc
from yt.utilities.exceptions import YTElementTypeNotRecognized
cdef extern from "mesh_triangulation.h":
@@ -54,6 +56,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef class LinearElementMesh:
@@ -240,39 +243,48 @@
'''
- cdef Patch* patches
cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef unsigned int mesh
# patches per element, vertices per element, vertices per face,
# and field points per element, respectively:
cdef int ppe, vpe, vpf, fpe
+ cdef int* faces
def __init__(self, YTEmbreeScene scene,
np.ndarray vertices,
np.ndarray indices,
np.ndarray field_data):
+ cdef int i, j
# 20-point hexes
if indices.shape[1] == 20:
self.vpe = 20
self.ppe = 6
self.vpf = 8
+ self.faces = <int*>malloc(6 * 8 * sizeof(int))
+ for i in range(6):
+ for j in range(8):
+ self.faces[i*6 + j] = hex20_faces[i][j]
+ self._build_from_indices_hex20(scene, vertices, indices, field_data)
# 10-point tets
elif indices.shape[1] == 10:
self.vpe = 10
self.ppe = 4
self.vpf = 6
+ self.faces = <int*>malloc(4 * 6 * sizeof(int))
+ for i in range(4):
+ for j in range(6):
+ self.faces[i*4 + j] = tet10_faces[i][j]
+ self._build_from_indices_tet10(scene, vertices, indices, field_data)
else:
raise NotImplementedError
- self._build_from_indices(scene, vertices, indices, field_data)
-
- cdef void _build_from_indices(self, YTEmbreeScene scene,
+ cdef void _build_from_indices_hex20(self, YTEmbreeScene scene,
np.ndarray vertices_in,
np.ndarray indices_in,
np.ndarray field_data):
- cdef int i, j, ind, idim
+ cdef int i, j, k, ind, idim
cdef int ne = indices_in.shape[0]
cdef int nv = vertices_in.shape[0]
cdef int npatch = self.ppe*ne;
@@ -283,14 +295,6 @@
self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
- cdef int faces[self.ppe][self.vpf];
- if self.vpe == 20:
- faces = hex20_faces
- elif self.vpe == 10:
- faces = tet10_faces
- else:
- raise NotImplementedError
-
for i in range(ne):
element_vertices = vertices_in[indices_in[i]]
for j in range(self.vpe):
@@ -305,22 +309,65 @@
patch = &(patches[i*self.ppe+j])
patch.geomID = mesh
for k in range(self.vpf): # for each vertex
- ind = faces[j][k]
+ ind = self.faces[j*self.ppe + k]
for idim in range(3): # for each spatial dimension (yikes)
patch.v[k][idim] = element_vertices[ind][idim]
patch.vertices = self.vertices + i*self.vpe*3
patch.field_data = self.field_data + i*self.vpe
- self.patches = patches
self.mesh = mesh
- rtcg.rtcSetUserData(scene.scene_i, self.mesh, self.patches)
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
<rtcgu.RTCBoundsFunc> patchBoundsFunc)
rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
<rtcgu.RTCIntersectFunc> patchIntersectFunc)
+ cdef void _build_from_indices_tet10(self, YTEmbreeScene scene,
+ np.ndarray vertices_in,
+ np.ndarray indices_in,
+ np.ndarray field_data):
+ cdef int i, j, k, ind, idim
+ cdef int ne = indices_in.shape[0]
+ cdef int nv = vertices_in.shape[0]
+ cdef int npatch = self.ppe*ne;
+
+ cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
+ cdef np.ndarray[np.float64_t, ndim=2] element_vertices
+ cdef Tet_Patch* tet_patches = <Tet_Patch*> malloc(npatch * sizeof(Tet_Patch))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
+
+ for i in range(ne):
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
+ for k in range(3):
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
+
+ cdef Tet_Patch* tet_patch
+ for i in range(ne): # for each element
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.ppe): # for each face
+ tet_patch = &(tet_patches[i*self.ppe+j])
+ tet_patch.geomID = mesh
+ for k in range(self.vpf): # for each vertex
+ ind = self.faces[j*self.ppe + k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ tet_patch.v[k][idim] = element_vertices[ind][idim]
+ tet_patch.vertices = self.vertices + i*self.vpe*3
+ tet_patch.field_data = self.field_data + i*self.vpe
+
+ self.mesh = mesh
+
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, tet_patches)
+ rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCBoundsFunc> tet_patchBoundsFunc)
+ rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCIntersectFunc> tet_patchIntersectFunc)
+
+
def __dealloc__(self):
- free(self.patches)
free(self.vertices)
free(self.field_data)
+ free(self.faces)
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_intersection.pxd
--- a/yt/utilities/lib/mesh_intersection.pxd
+++ b/yt/utilities/lib/mesh_intersection.pxd
@@ -1,7 +1,7 @@
cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry as rtcg
-from yt.utilities.lib.mesh_construction cimport Patch
+from yt.utilities.lib.mesh_construction cimport Patch, Tet_Patch
cimport cython
cdef void patchIntersectFunc(Patch* patches,
@@ -11,3 +11,11 @@
cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil
+
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil
+
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_intersection.pyx
--- a/yt/utilities/lib/mesh_intersection.pyx
+++ b/yt/utilities/lib/mesh_intersection.pyx
@@ -20,26 +20,30 @@
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
-from yt.utilities.lib.mesh_samplers cimport sample_hex20
+from yt.utilities.lib.mesh_samplers cimport sample_hex20, sample_tet10
from yt.utilities.lib.bounding_volume_hierarchy cimport BBox
from yt.utilities.lib.primitives cimport \
patchSurfaceFunc, \
patchSurfaceDerivU, \
patchSurfaceDerivV, \
RayHitData, \
- compute_patch_hit
+ compute_patch_hit, \
+ tet_patchSurfaceFunc, \
+ tet_patchSurfaceDerivU, \
+ tet_patchSurfaceDerivV, \
+ compute_tet_patch_hit
from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void patchBoundsFunc(Patch* patches,
+cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:
cdef Patch patch = patches[item]
-
+
cdef float lo_x = 1.0e300
cdef float lo_y = 1.0e300
cdef float lo_z = 1.0e300
@@ -88,8 +92,71 @@
ray.geomID = patch.geomID
ray.primID = item
ray.Ng[0] = hd.t
-
+
# sample the solution at the calculated point
sample_hex20(patches, ray)
return
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef float lo_x = 1.0e300
+ cdef float lo_y = 1.0e300
+ cdef float lo_z = 1.0e300
+
+ cdef float hi_x = -1.0e300
+ cdef float hi_y = -1.0e300
+ cdef float hi_z = -1.0e300
+
+ cdef int i
+ for i in range(6):
+ lo_x = fmin(lo_x, tet_patch.v[i][0])
+ lo_y = fmin(lo_y, tet_patch.v[i][1])
+ lo_z = fmin(lo_z, tet_patch.v[i][2])
+ hi_x = fmax(hi_x, tet_patch.v[i][0])
+ hi_y = fmax(hi_y, tet_patch.v[i][1])
+ hi_z = fmax(hi_z, tet_patch.v[i][2])
+
+ bounds_o.lower_x = lo_x
+ bounds_o.lower_y = lo_y
+ bounds_o.lower_z = lo_z
+ bounds_o.upper_x = hi_x
+ bounds_o.upper_y = hi_y
+ bounds_o.upper_z = hi_z
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.org, ray.dir)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):
+ return
+
+ if (hd.converged and 0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1):
+
+ # we have a hit, so update ray information
+ ray.u = hd.u
+ ray.v = hd.v
+ ray.geomID = tet_patch.geomID
+ ray.primID = item
+ ray.Ng[0] = hd.t
+
+ # sample the solution at the calculated point
+ sample_tet10(tet_patches, ray)
+
+ return
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_samplers.pxd
--- a/yt/utilities/lib/mesh_samplers.pxd
+++ b/yt/utilities/lib/mesh_samplers.pxd
@@ -13,3 +13,6 @@
cdef void sample_hex20(void* userPtr,
rtcr.RTCRay& ray) nogil
+
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil
diff -r 897ffcfa02ef845f1f3dec541893e0e45e4cdefe -r f42ee7dc1a8879a528678b724735484d1e51a822 yt/utilities/lib/mesh_samplers.pyx
--- a/yt/utilities/lib/mesh_samplers.pyx
+++ b/yt/utilities/lib/mesh_samplers.pyx
@@ -18,14 +18,16 @@
from pyembree.rtcore cimport Vec3f, Triangle, Vertex
from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
- Patch
-from yt.utilities.lib.primitives cimport patchSurfaceFunc
+ Patch, \
+ Tet_Patch
+from yt.utilities.lib.primitives cimport patchSurfaceFunc, tet_patchSurfaceFunc
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \
Q1Sampler3D, \
S2Sampler3D, \
- W1Sampler3D
+ W1Sampler3D, \
+ Tet2Sampler3D
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmax
@@ -34,6 +36,7 @@
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
@cython.boundscheck(False)
@@ -94,7 +97,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(8):
element_indices[i] = data.element_indices[elem_id*8+i]
@@ -104,7 +107,7 @@
for i in range(8):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -145,7 +148,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(6):
element_indices[i] = data.element_indices[elem_id*6+i]
@@ -155,7 +158,7 @@
for i in range(6):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -195,14 +198,14 @@
patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
for i in range(3):
position[i] = <double> pos[i]
-
+
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
S2Sampler.map_real_to_unit(mapped_coord, patch.vertices, position)
val = S2Sampler.sample_at_unit_point(mapped_coord, patch.field_data)
ray.time = val
ray.instID = S2Sampler.check_mesh_lines(mapped_coord)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -228,14 +231,14 @@
# ray_id records the id number of the hit according to
# embree, in which the primitives are triangles. Here,
# we convert this to the element id by dividing by the
- # number of triangles per element.
+ # number of triangles per element.
elem_id = ray_id / data.tpe
for i in range(4):
element_indices[i] = data.element_indices[elem_id*4+i]
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
for i in range(data.fpe):
field_data[i] = data.field_data[elem_id*data.fpe+i]
@@ -249,3 +252,39 @@
val = P1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
ray.instID = P1Sampler.check_mesh_lines(mapped_coord)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.initializedcheck(False)
+ at cython.cdivision(True)
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil:
+ cdef int ray_id, elem_id, i
+ cdef double val
+ cdef double[3] position
+ cdef float[3] pos
+ cdef Tet_Patch* data
+
+ data = <Tet_Patch*> userPtr
+ ray_id = ray.primID
+ if ray_id == -1:
+ return
+ cdef Tet_Patch tet_patch = data[ray_id]
+
+ # ray_id records the id number of the hit according to
+ # embree, in which the primitives are patches. Here,
+ # we convert this to the element id by dividing by the
+ # number of patches per element.
+ elem_id = ray_id / 4
+
+ # fills "position" with the physical position of the hit
+ tet_patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
+ for i in range(3):
+ position[i] = <double> pos[i]
+
+ # we use ray.time to pass the value of the field
+ cdef double mapped_coord[3]
+ Tet2Sampler.map_real_to_unit(mapped_coord, tet_patch.vertices, position)
+ val = Tet2Sampler.sample_at_unit_point(mapped_coord, tet_patch.field_data)
+ ray.time = val
+ ray.instID = Tet2Sampler.check_mesh_lines(mapped_coord)
https://bitbucket.org/yt_analysis/yt/commits/4783f44f889e/
Changeset: 4783f44f889e
Branch: yt
User: al007
Date: 2016-10-05 21:43:14+00:00
Summary: Fixed embree rendering.
Affected #: 3 files
diff -r f42ee7dc1a8879a528678b724735484d1e51a822 -r 4783f44f889e249cde5ae0c7ef2cfd823216d292 setup.py
--- a/setup.py
+++ b/setup.py
@@ -322,7 +322,7 @@
def finalize_options(self):
from Cython.Build import cythonize
self.distribution.ext_modules[:] = cythonize(
- self.distribution.ext_modules, gdb_debug=True)
+ self.distribution.ext_modules)
_build_ext.finalize_options(self)
# Prevent numpy from thinking it is still in its setup process
# see http://stackoverflow.com/a/21621493/1382869
@@ -337,7 +337,7 @@
def run(self):
# Make sure the compiled Cython files in the distribution are up-to-date
from Cython.Build import cythonize
- cythonize(cython_extensions, gdb_dbg=True)
+ cythonize(cython_extensions)
_sdist.run(self)
setup(
diff -r f42ee7dc1a8879a528678b724735484d1e51a822 -r 4783f44f889e249cde5ae0c7ef2cfd823216d292 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -900,6 +900,31 @@
return 0
return 1
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ cdef double u, v
+ cdef double thresh = 2.0e-2
+ if mapped_coord[0] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ elif mapped_coord[1] == 0:
+ u = mapped_coord[2]
+ v = mapped_coord[0]
+ elif mapped_coord[2] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[0]
+ else:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ if ((u < thresh) or
+ (v < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh)):
+ return 1
+ return -1
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r f42ee7dc1a8879a528678b724735484d1e51a822 -r 4783f44f889e249cde5ae0c7ef2cfd823216d292 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -249,7 +249,6 @@
# patches per element, vertices per element, vertices per face,
# and field points per element, respectively:
cdef int ppe, vpe, vpf, fpe
- cdef int* faces
def __init__(self, YTEmbreeScene scene,
np.ndarray vertices,
@@ -262,20 +261,12 @@
self.vpe = 20
self.ppe = 6
self.vpf = 8
- self.faces = <int*>malloc(6 * 8 * sizeof(int))
- for i in range(6):
- for j in range(8):
- self.faces[i*6 + j] = hex20_faces[i][j]
self._build_from_indices_hex20(scene, vertices, indices, field_data)
# 10-point tets
elif indices.shape[1] == 10:
self.vpe = 10
self.ppe = 4
self.vpf = 6
- self.faces = <int*>malloc(4 * 6 * sizeof(int))
- for i in range(4):
- for j in range(6):
- self.faces[i*4 + j] = tet10_faces[i][j]
self._build_from_indices_tet10(scene, vertices, indices, field_data)
else:
raise NotImplementedError
@@ -309,7 +300,7 @@
patch = &(patches[i*self.ppe+j])
patch.geomID = mesh
for k in range(self.vpf): # for each vertex
- ind = self.faces[j*self.ppe + k]
+ ind = hex20_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
patch.v[k][idim] = element_vertices[ind][idim]
patch.vertices = self.vertices + i*self.vpe*3
@@ -352,7 +343,7 @@
tet_patch = &(tet_patches[i*self.ppe+j])
tet_patch.geomID = mesh
for k in range(self.vpf): # for each vertex
- ind = self.faces[j*self.ppe + k]
+ ind = tet10_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
tet_patch.v[k][idim] = element_vertices[ind][idim]
tet_patch.vertices = self.vertices + i*self.vpe*3
@@ -370,4 +361,3 @@
def __dealloc__(self):
free(self.vertices)
free(self.field_data)
- free(self.faces)
https://bitbucket.org/yt_analysis/yt/commits/62f20c5cdc2b/
Changeset: 62f20c5cdc2b
Branch: yt
User: al007
Date: 2016-10-06 17:23:06+00:00
Summary: Add second order tetrahedral answer test for volume rendering.
Affected #: 1 file
diff -r 4783f44f889e249cde5ae0c7ef2cfd823216d292 -r 62f20c5cdc2b881d197b8e3ca0d1669f28921f8e yt/visualization/volume_rendering/tests/test_mesh_render.py
--- a/yt/visualization/volume_rendering/tests/test_mesh_render.py
+++ b/yt/visualization/volume_rendering/tests/test_mesh_render.py
@@ -55,7 +55,7 @@
images.append(im)
return images
-
+
@requires_module("pyembree")
def test_surface_mesh_render_pyembree():
ytcfg["yt", "ray_tracing_engine"] = "embree"
@@ -145,7 +145,7 @@
im = sc.render()
return compare(ds, im, "%s_render_answers_wedge6_%s_%s" %
(engine, field[0], field[1]))
-
+
@requires_ds(wedge6)
@requires_module("pyembree")
def test_wedge6_render_pyembree():
@@ -157,6 +157,28 @@
for field in wedge6_fields:
yield wedge6_render("yt", field)
+tet10 = "SecondOrderTets/tet10_unstructured_out.e"
+tet10_fields = [('connect1', 'uz')]
+
+def tet10_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(tet10, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_tet10_%s_%s" %
+ (engine, field[0], field[1]))
+
+ at requires_ds(tet10)
+ at requires_module("pyembree")
+def test_tet10_render_pyembree():
+ for field in tet10_fields:
+ yield tet10_render("embree", field)
+
+ at requires_ds(tet10)
+def test_tet10_render():
+ for field in tet10_fields:
+ yield tet10_render("yt", field)
+
def perspective_mesh_render(engine):
ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(hex8)
@@ -169,7 +191,7 @@
cam.resolution = (800, 800)
im = sc.render()
return compare(ds, im, "%s_perspective_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_perspective_mesh_render_pyembree():
@@ -195,7 +217,7 @@
sc.add_source(ms2)
im = sc.render()
return compare(ds, im, "%s_composite_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_composite_mesh_render_pyembree():
https://bitbucket.org/yt_analysis/yt/commits/ab97afcce237/
Changeset: ab97afcce237
Branch: yt
User: al007
Date: 2016-10-06 19:05:19+00:00
Summary: Add void* patches class member so we can dealloc later. Increment test.
Affected #: 2 files
diff -r 62f20c5cdc2b881d197b8e3ca0d1669f28921f8e -r ab97afcce237b906bd308f6ad6961f58ca6b1ddd tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -50,7 +50,7 @@
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_005:
+ local_varia_006:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
diff -r 62f20c5cdc2b881d197b8e3ca0d1669f28921f8e -r ab97afcce237b906bd308f6ad6961f58ca6b1ddd yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -243,6 +243,7 @@
'''
+ cdef void* patches
cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef unsigned int mesh
@@ -325,7 +326,7 @@
cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
cdef np.ndarray[np.float64_t, ndim=2] element_vertices
- cdef Tet_Patch* tet_patches = <Tet_Patch*> malloc(npatch * sizeof(Tet_Patch))
+ cdef Tet_Patch* patches = <Tet_Patch*> malloc(npatch * sizeof(Tet_Patch))
self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
@@ -336,22 +337,22 @@
for k in range(3):
self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
- cdef Tet_Patch* tet_patch
+ cdef Tet_Patch* patch
for i in range(ne): # for each element
element_vertices = vertices_in[indices_in[i]]
for j in range(self.ppe): # for each face
- tet_patch = &(tet_patches[i*self.ppe+j])
- tet_patch.geomID = mesh
+ patch = &(patches[i*self.ppe+j])
+ patch.geomID = mesh
for k in range(self.vpf): # for each vertex
ind = tet10_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
- tet_patch.v[k][idim] = element_vertices[ind][idim]
- tet_patch.vertices = self.vertices + i*self.vpe*3
- tet_patch.field_data = self.field_data + i*self.vpe
+ patch.v[k][idim] = element_vertices[ind][idim]
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
self.mesh = mesh
- rtcg.rtcSetUserData(scene.scene_i, self.mesh, tet_patches)
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
<rtcgu.RTCBoundsFunc> tet_patchBoundsFunc)
rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
@@ -359,5 +360,6 @@
def __dealloc__(self):
+ free(self.patches)
free(self.vertices)
free(self.field_data)
https://bitbucket.org/yt_analysis/yt/commits/c84972e66800/
Changeset: c84972e66800
Branch: yt
User: al007
Date: 2016-10-07 14:03:52+00:00
Summary: Assign patches in _build_indices to self.patches for memory deallocation.
Affected #: 1 file
diff -r ab97afcce237b906bd308f6ad6961f58ca6b1ddd -r c84972e66800d18539998c65871828e25e6fb344 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -307,6 +307,7 @@
patch.vertices = self.vertices + i*self.vpe*3
patch.field_data = self.field_data + i*self.vpe
+ self.patches = patches
self.mesh = mesh
rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
@@ -350,6 +351,7 @@
patch.vertices = self.vertices + i*self.vpe*3
patch.field_data = self.field_data + i*self.vpe
+ self.patches = patches
self.mesh = mesh
rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
https://bitbucket.org/yt_analysis/yt/commits/1d43d35f96cf/
Changeset: 1d43d35f96cf
Branch: yt
User: ngoldbaum
Date: 2016-10-17 21:19:44+00:00
Summary: Merged in al007/yt (pull request #2401)
Implement ability to volume render 3D meshes composed of second-order tetrahedrons
Affected #: 20 files
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -50,7 +50,7 @@
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_005:
+ local_varia_006:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
# node for the bounding volume hierarchy
cdef struct BVHNode:
@@ -69,8 +70,11 @@
cdef void _set_up_patches(self,
np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil
+ cdef void _set_up_tet_patches(self,
+ np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil
cdef void intersect(self, Ray* ray) nogil
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil
cdef void _recursive_intersect(self, Ray* ray, BVHNode* node) nogil
cdef BVHNode* _recursive_build(self, np.int64_t begin, np.int64_t end) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -18,19 +18,27 @@
Patch, \
ray_patch_intersect, \
patch_centroid, \
- patch_bbox
+ patch_bbox, \
+ TetPatch, \
+ ray_tet_patch_intersect, \
+ tet_patch_centroid, \
+ tet_patch_bbox
+
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
P1Sampler3D, \
W1Sampler3D, \
- S2Sampler3D
+ S2Sampler3D, \
+ Tet2Sampler3D
+
from yt.utilities.lib.vec3_ops cimport L2_norm
cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
@@ -45,11 +53,11 @@
'''
This class implements a bounding volume hierarchy (BVH), a spatial acceleration
- structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
- partitioning the *volume* of the parent to create the children, we partition the
+ structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
+ partitioning the *volume* of the parent to create the children, we partition the
triangles themselves into 'left' or 'right' sub-trees. The bounding volume for a
node is then determined by computing the bounding volume of the triangles that
- belong to it. This allows us to quickly discard triangles that are not close
+ belong to it. This allows us to quickly discard triangles that are not close
to intersecting a given ray.
'''
@@ -82,6 +90,9 @@
elif self.num_verts_per_elem == 20:
self.num_prim_per_elem = 6
self.sampler = S2Sampler
+ elif self.num_verts_per_elem == 10:
+ self.num_prim_per_elem = 4
+ self.sampler = Tet2Sampler
else:
raise NotImplementedError("Could not determine element type for "
"nverts = %d. " % self.num_verts_per_elem)
@@ -109,7 +120,7 @@
self.vertices[vertex_offset + k] = vertices[indices[i,j]][k]
field_offset = i*self.num_field_per_elem
for j in range(self.num_field_per_elem):
- self.field_data[field_offset + j] = field_data[i][j]
+ self.field_data[field_offset + j] = field_data[i][j]
# set up primitives
if self.num_verts_per_elem == 20:
@@ -118,13 +129,19 @@
self.get_bbox = patch_bbox
self.get_intersect = ray_patch_intersect
self._set_up_patches(vertices, indices)
+ elif self.num_verts_per_elem == 10:
+ self.primitives = malloc(self.num_prim * sizeof(TetPatch))
+ self.get_centroid = tet_patch_centroid
+ self.get_bbox = tet_patch_bbox
+ self.get_intersect = ray_tet_patch_intersect
+ self._set_up_tet_patches(vertices, indices)
else:
self.primitives = malloc(self.num_prim * sizeof(Triangle))
self.get_centroid = triangle_centroid
self.get_bbox = triangle_bbox
self.get_intersect = ray_triangle_intersect
self._set_up_triangles(vertices, indices)
-
+
self.root = self._recursive_build(0, self.num_prim)
@cython.boundscheck(False)
@@ -156,6 +173,32 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+ cdef void _set_up_tet_patches(self, np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil:
+ cdef TetPatch* tet_patch
+ cdef np.int64_t i, j, k, ind, idim
+ cdef np.int64_t offset, prim_index
+ for i in range(self.num_elem):
+ offset = self.num_prim_per_elem*i
+ for j in range(self.num_prim_per_elem): # for each face
+ prim_index = offset + j
+ tet_patch = &( <TetPatch*> self.primitives)[prim_index]
+ self.prim_ids[prim_index] = prim_index
+ tet_patch.elem_id = i
+ for k in range(6): # for each vertex
+ ind = tet10_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ tet_patch.v[k][idim] = vertices[indices[i, ind]][idim]
+ self.get_centroid(self.primitives,
+ prim_index,
+ self.centroids[prim_index])
+ self.get_bbox(self.primitives,
+ prim_index,
+ &(self.bboxes[prim_index]))
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
cdef void _set_up_triangles(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:
# fill our array of primitives
@@ -181,7 +224,7 @@
tri_index,
self.centroids[tri_index])
self.get_bbox(self.primitives,
- tri_index,
+ tri_index,
&(self.bboxes[tri_index]))
cdef void _recursive_free(self, BVHNode* node) nogil:
@@ -191,6 +234,8 @@
free(node)
def __dealloc__(self):
+ if self.root == NULL:
+ return
self._recursive_free(self.root)
free(self.primitives)
free(self.prim_ids)
@@ -224,11 +269,11 @@
mid += 1
begin += 1
return mid
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil:
cdef np.int64_t i, j
cdef BBox box = self.bboxes[begin]
@@ -236,7 +281,7 @@
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],
self.bboxes[i].left_edge[j])
- box.right_edge[j] = fmax(box.right_edge[j],
+ box.right_edge[j] = fmax(box.right_edge[j],
self.bboxes[i].right_edge[j])
node.bbox = box
@@ -245,7 +290,7 @@
@cython.cdivision(True)
cdef void intersect(self, Ray* ray) nogil:
self._recursive_intersect(ray, self.root)
-
+
if ray.elem_id < 0:
return
@@ -253,7 +298,7 @@
cdef np.int64_t i
for i in range(3):
position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
-
+
cdef np.float64_t* vertex_ptr
cdef np.float64_t* field_ptr
vertex_ptr = self.vertices + ray.elem_id*self.num_verts_per_elem*3
@@ -297,17 +342,17 @@
node.end = end
self._get_node_bbox(node, begin, end)
-
+
# check for leaf
if (end - begin) <= LEAF_SIZE:
return node
-
+
# we use the "split in the middle of the longest axis approach"
# see: http://www.vadimkravcenko.com/bvh-tree-building/
# compute longest dimension
cdef np.int64_t ax = 0
- cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
+ cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
node.bbox.left_edge[0])
if fabs(node.bbox.right_edge[1] - node.bbox.left_edge[1]) > d:
ax = 1
@@ -323,7 +368,7 @@
if(mid == begin or mid == end):
mid = begin + (end-begin)/2
-
+
# recursively build sub-trees
node.left = self._recursive_build(begin, mid)
node.right = self._recursive_build(mid, end)
@@ -337,16 +382,16 @@
cdef void cast_rays(np.float64_t* image,
const np.float64_t* origins,
const np.float64_t* direction,
- const int N,
+ const int N,
BVH bvh) nogil:
- cdef Ray* ray
+ cdef Ray* ray
cdef int i, j, k
-
+
with nogil, parallel():
-
+
ray = <Ray *> malloc(sizeof(Ray))
-
+
for k in range(3):
ray.direction[k] = direction[k]
ray.inv_dir[k] = 1.0 / direction[k]
@@ -366,11 +411,11 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
+def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
np.ndarray[np.float64_t, ndim=2] origins,
np.ndarray[np.float64_t, ndim=1] direction,
BVH bvh):
-
+
cdef int N = origins.shape[0]
cast_rays(&image[0], &origins[0, 0], &direction[0], N, bvh)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -834,11 +834,10 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tris, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(2):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] - 1.0 > self.inclusion_tol):
+ return 0
return 1
cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
@@ -893,13 +892,39 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tets, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(3):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[2] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] + mapped_coord[2] - 1.0 > \
+ self.inclusion_tol):
+ return 0
return 1
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ cdef double u, v
+ cdef double thresh = 2.0e-2
+ if mapped_coord[0] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ elif mapped_coord[1] == 0:
+ u = mapped_coord[2]
+ v = mapped_coord[0]
+ elif mapped_coord[2] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[0]
+ else:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ if ((u < thresh) or
+ (v < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh)):
+ return 1
+ return -1
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_construction.pxd
--- a/yt/utilities/lib/mesh_construction.pxd
+++ b/yt/utilities/lib/mesh_construction.pxd
@@ -18,3 +18,9 @@
unsigned int geomID
np.float64_t* vertices
np.float64_t* field_data
+
+ctypedef struct Tet_Patch:
+ float[6][3] v
+ unsigned int geomID
+ np.float64_t* vertices
+ np.float64_t* field_data
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -1,5 +1,5 @@
"""
-This file contains the ElementMesh classes, which represent the target that the
+This file contains the ElementMesh classes, which represent the target that the
rays will be cast at when rendering finite element data. This class handles
the interface between the internal representation of the mesh and the pyembree
representation.
@@ -21,7 +21,7 @@
from libc.math cimport fmax, sqrt
cimport numpy as np
-cimport pyembree.rtcore as rtc
+cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu
@@ -37,13 +37,15 @@
sample_wedge
from mesh_intersection cimport \
patchIntersectFunc, \
- patchBoundsFunc
+ patchBoundsFunc, \
+ tet_patchIntersectFunc, \
+ tet_patchBoundsFunc
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
@@ -54,13 +56,14 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef class LinearElementMesh:
r'''
This creates a 1st-order mesh to be ray-traced with embree.
- Currently, we handle non-triangular mesh types by converting them
+ Currently, we handle non-triangular mesh types by converting them
to triangular meshes. This class performs this transformation.
Currently, this is implemented for hexahedral and tetrahedral
meshes.
@@ -71,21 +74,21 @@
scene : EmbreeScene
This is the scene to which the constructed polygons will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the polygon mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the polygon mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
- This should either have the shape (num_elements, 4) or
- (num_elements, 8) for tetrahedral and hexahedral meshes,
- respectively. For tetrahedral meshes, each element will
+ This should either have the shape (num_elements, 4) or
+ (num_elements, 8) for tetrahedral and hexahedral meshes,
+ respectively. For tetrahedral meshes, each element will
be represented by four triangles in the scene. For hex meshes,
- each element will be represented by 12 triangles, 2 for each
+ each element will be represented by 12 triangles, 2 for each
face. For hex meshes, we assume that the node ordering is as
- defined here:
+ defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
cdef Vertex* vertices
@@ -93,7 +96,7 @@
cdef unsigned int mesh
cdef double* field_data
cdef rtcg.RTCFilterFunc filter_func
- # triangles per element, vertices per element, and field points per
+ # triangles per element, vertices per element, and field points per
# element, respectively:
cdef int tpe, vpe, fpe
cdef int[MAX_NUM_TRI][3] tri_array
@@ -101,7 +104,7 @@
cdef MeshDataContainer datac
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray data):
@@ -119,7 +122,7 @@
self.tpe = TETRA_NT
self.tri_array = triangulate_tetra
else:
- raise YTElementTypeNotRecognized(vertices.shape[1],
+ raise YTElementTypeNotRecognized(vertices.shape[1],
indices.shape[1])
self._build_from_indices(scene, vertices, indices)
@@ -142,7 +145,7 @@
for i in range(nv):
vertices[i].x = vertices_in[i, 0]
vertices[i].y = vertices_in[i, 1]
- vertices[i].z = vertices_in[i, 2]
+ vertices[i].z = vertices_in[i, 2]
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_VERTEX_BUFFER,
vertices, 0, sizeof(Vertex))
@@ -156,7 +159,7 @@
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_INDEX_BUFFER,
triangles, 0, sizeof(Triangle))
- cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
+ cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
for i in range(ne):
for j in range(self.vpe):
element_indices[i*self.vpe + j] = indices_in[i][j]
@@ -189,7 +192,7 @@
datac.vpe = self.vpe
datac.fpe = self.fpe
self.datac = datac
-
+
rtcg.rtcSetUserData(scene.scene_i, self.mesh, &self.datac)
cdef void _set_sampler_type(self, YTEmbreeScene scene):
@@ -206,7 +209,7 @@
rtcg.rtcSetIntersectionFilterFunction(scene.scene_i,
self.mesh,
self.filter_func)
-
+
def __dealloc__(self):
free(self.field_data)
free(self.element_indices)
@@ -227,84 +230,137 @@
scene : EmbreeScene
This is the scene to which the constructed patches will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
This should have the shape (num_elements, 20). Each hex will be
- represented in the scene by 6 bi-quadratic patches. We assume that
- the node ordering is as defined here:
+ represented in the scene by 6 bi-quadratic patches. We assume that
+ the node ordering is as defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
- cdef Patch* patches
+ cdef void* patches
cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef unsigned int mesh
- # patches per element, vertices per element, and field points per
- # element, respectively:
- cdef int ppe, vpe, fpe
+ # patches per element, vertices per element, vertices per face,
+ # and field points per element, respectively:
+ cdef int ppe, vpe, vpf, fpe
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray field_data):
+ cdef int i, j
- # only 20-point hexes are supported right now.
+ # 20-point hexes
if indices.shape[1] == 20:
self.vpe = 20
+ self.ppe = 6
+ self.vpf = 8
+ self._build_from_indices_hex20(scene, vertices, indices, field_data)
+ # 10-point tets
+ elif indices.shape[1] == 10:
+ self.vpe = 10
+ self.ppe = 4
+ self.vpf = 6
+ self._build_from_indices_tet10(scene, vertices, indices, field_data)
else:
raise NotImplementedError
- self._build_from_indices(scene, vertices, indices, field_data)
-
- cdef void _build_from_indices(self, YTEmbreeScene scene,
+ cdef void _build_from_indices_hex20(self, YTEmbreeScene scene,
np.ndarray vertices_in,
np.ndarray indices_in,
np.ndarray field_data):
- cdef int i, j, ind, idim
+ cdef int i, j, k, ind, idim
cdef int ne = indices_in.shape[0]
cdef int nv = vertices_in.shape[0]
- cdef int npatch = 6*ne;
+ cdef int npatch = self.ppe*ne;
cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
cdef np.ndarray[np.float64_t, ndim=2] element_vertices
cdef Patch* patches = <Patch*> malloc(npatch * sizeof(Patch))
- self.vertices = <np.float64_t*> malloc(20 * ne * 3 * sizeof(np.float64_t))
- self.field_data = <np.float64_t*> malloc(20 * ne * sizeof(np.float64_t))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
for i in range(ne):
element_vertices = vertices_in[indices_in[i]]
- for j in range(20):
- self.field_data[i*20 + j] = field_data[i][j]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
for k in range(3):
- self.vertices[i*20*3 + j*3 + k] = element_vertices[j][k]
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
cdef Patch* patch
for i in range(ne): # for each element
element_vertices = vertices_in[indices_in[i]]
- for j in range(6): # for each face
- patch = &(patches[i*6+j])
+ for j in range(self.ppe): # for each face
+ patch = &(patches[i*self.ppe+j])
patch.geomID = mesh
- for k in range(8): # for each vertex
+ for k in range(self.vpf): # for each vertex
ind = hex20_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
patch.v[k][idim] = element_vertices[ind][idim]
- patch.vertices = self.vertices + i*20*3
- patch.field_data = self.field_data + i*20
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
self.patches = patches
self.mesh = mesh
- rtcg.rtcSetUserData(scene.scene_i, self.mesh, self.patches)
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
<rtcgu.RTCBoundsFunc> patchBoundsFunc)
rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
<rtcgu.RTCIntersectFunc> patchIntersectFunc)
+ cdef void _build_from_indices_tet10(self, YTEmbreeScene scene,
+ np.ndarray vertices_in,
+ np.ndarray indices_in,
+ np.ndarray field_data):
+ cdef int i, j, k, ind, idim
+ cdef int ne = indices_in.shape[0]
+ cdef int nv = vertices_in.shape[0]
+ cdef int npatch = self.ppe*ne;
+
+ cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
+ cdef np.ndarray[np.float64_t, ndim=2] element_vertices
+ cdef Tet_Patch* patches = <Tet_Patch*> malloc(npatch * sizeof(Tet_Patch))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
+
+ for i in range(ne):
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
+ for k in range(3):
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
+
+ cdef Tet_Patch* patch
+ for i in range(ne): # for each element
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.ppe): # for each face
+ patch = &(patches[i*self.ppe+j])
+ patch.geomID = mesh
+ for k in range(self.vpf): # for each vertex
+ ind = tet10_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ patch.v[k][idim] = element_vertices[ind][idim]
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
+
+ self.patches = patches
+ self.mesh = mesh
+
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
+ rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCBoundsFunc> tet_patchBoundsFunc)
+ rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCIntersectFunc> tet_patchIntersectFunc)
+
+
def __dealloc__(self):
free(self.patches)
free(self.vertices)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_intersection.pxd
--- a/yt/utilities/lib/mesh_intersection.pxd
+++ b/yt/utilities/lib/mesh_intersection.pxd
@@ -1,7 +1,7 @@
cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry as rtcg
-from yt.utilities.lib.mesh_construction cimport Patch
+from yt.utilities.lib.mesh_construction cimport Patch, Tet_Patch
cimport cython
cdef void patchIntersectFunc(Patch* patches,
@@ -11,3 +11,11 @@
cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil
+
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil
+
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_intersection.pyx
--- a/yt/utilities/lib/mesh_intersection.pyx
+++ b/yt/utilities/lib/mesh_intersection.pyx
@@ -20,26 +20,30 @@
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
-from yt.utilities.lib.mesh_samplers cimport sample_hex20
+from yt.utilities.lib.mesh_samplers cimport sample_hex20, sample_tet10
from yt.utilities.lib.bounding_volume_hierarchy cimport BBox
from yt.utilities.lib.primitives cimport \
patchSurfaceFunc, \
patchSurfaceDerivU, \
patchSurfaceDerivV, \
RayHitData, \
- compute_patch_hit
+ compute_patch_hit, \
+ tet_patchSurfaceFunc, \
+ tet_patchSurfaceDerivU, \
+ tet_patchSurfaceDerivV, \
+ compute_tet_patch_hit
from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void patchBoundsFunc(Patch* patches,
+cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:
cdef Patch patch = patches[item]
-
+
cdef float lo_x = 1.0e300
cdef float lo_y = 1.0e300
cdef float lo_z = 1.0e300
@@ -88,8 +92,71 @@
ray.geomID = patch.geomID
ray.primID = item
ray.Ng[0] = hd.t
-
+
# sample the solution at the calculated point
sample_hex20(patches, ray)
return
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef float lo_x = 1.0e300
+ cdef float lo_y = 1.0e300
+ cdef float lo_z = 1.0e300
+
+ cdef float hi_x = -1.0e300
+ cdef float hi_y = -1.0e300
+ cdef float hi_z = -1.0e300
+
+ cdef int i
+ for i in range(6):
+ lo_x = fmin(lo_x, tet_patch.v[i][0])
+ lo_y = fmin(lo_y, tet_patch.v[i][1])
+ lo_z = fmin(lo_z, tet_patch.v[i][2])
+ hi_x = fmax(hi_x, tet_patch.v[i][0])
+ hi_y = fmax(hi_y, tet_patch.v[i][1])
+ hi_z = fmax(hi_z, tet_patch.v[i][2])
+
+ bounds_o.lower_x = lo_x
+ bounds_o.lower_y = lo_y
+ bounds_o.lower_z = lo_z
+ bounds_o.upper_x = hi_x
+ bounds_o.upper_y = hi_y
+ bounds_o.upper_z = hi_z
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.org, ray.dir)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):
+ return
+
+ if (hd.converged and 0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1):
+
+ # we have a hit, so update ray information
+ ray.u = hd.u
+ ray.v = hd.v
+ ray.geomID = tet_patch.geomID
+ ray.primID = item
+ ray.Ng[0] = hd.t
+
+ # sample the solution at the calculated point
+ sample_tet10(tet_patches, ray)
+
+ return
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_samplers.pxd
--- a/yt/utilities/lib/mesh_samplers.pxd
+++ b/yt/utilities/lib/mesh_samplers.pxd
@@ -13,3 +13,6 @@
cdef void sample_hex20(void* userPtr,
rtcr.RTCRay& ray) nogil
+
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_samplers.pyx
--- a/yt/utilities/lib/mesh_samplers.pyx
+++ b/yt/utilities/lib/mesh_samplers.pyx
@@ -18,14 +18,16 @@
from pyembree.rtcore cimport Vec3f, Triangle, Vertex
from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
- Patch
-from yt.utilities.lib.primitives cimport patchSurfaceFunc
+ Patch, \
+ Tet_Patch
+from yt.utilities.lib.primitives cimport patchSurfaceFunc, tet_patchSurfaceFunc
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \
Q1Sampler3D, \
S2Sampler3D, \
- W1Sampler3D
+ W1Sampler3D, \
+ Tet2Sampler3D
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmax
@@ -34,6 +36,7 @@
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
@cython.boundscheck(False)
@@ -94,7 +97,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(8):
element_indices[i] = data.element_indices[elem_id*8+i]
@@ -104,7 +107,7 @@
for i in range(8):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -145,7 +148,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(6):
element_indices[i] = data.element_indices[elem_id*6+i]
@@ -155,7 +158,7 @@
for i in range(6):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -195,14 +198,14 @@
patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
for i in range(3):
position[i] = <double> pos[i]
-
+
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
S2Sampler.map_real_to_unit(mapped_coord, patch.vertices, position)
val = S2Sampler.sample_at_unit_point(mapped_coord, patch.field_data)
ray.time = val
ray.instID = S2Sampler.check_mesh_lines(mapped_coord)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -228,14 +231,14 @@
# ray_id records the id number of the hit according to
# embree, in which the primitives are triangles. Here,
# we convert this to the element id by dividing by the
- # number of triangles per element.
+ # number of triangles per element.
elem_id = ray_id / data.tpe
for i in range(4):
element_indices[i] = data.element_indices[elem_id*4+i]
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
for i in range(data.fpe):
field_data[i] = data.field_data[elem_id*data.fpe+i]
@@ -249,3 +252,39 @@
val = P1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
ray.instID = P1Sampler.check_mesh_lines(mapped_coord)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.initializedcheck(False)
+ at cython.cdivision(True)
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil:
+ cdef int ray_id, elem_id, i
+ cdef double val
+ cdef double[3] position
+ cdef float[3] pos
+ cdef Tet_Patch* data
+
+ data = <Tet_Patch*> userPtr
+ ray_id = ray.primID
+ if ray_id == -1:
+ return
+ cdef Tet_Patch tet_patch = data[ray_id]
+
+ # ray_id records the id number of the hit according to
+ # embree, in which the primitives are patches. Here,
+ # we convert this to the element id by dividing by the
+ # number of patches per element.
+ elem_id = ray_id / 4
+
+ # fills "position" with the physical position of the hit
+ tet_patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
+ for i in range(3):
+ position[i] = <double> pos[i]
+
+ # we use ray.time to pass the value of the field
+ cdef double mapped_coord[3]
+ Tet2Sampler.map_real_to_unit(mapped_coord, tet_patch.vertices, position)
+ val = Tet2Sampler.sample_at_unit_point(mapped_coord, tet_patch.field_data)
+ ray.time = val
+ ray.instID = Tet2Sampler.check_mesh_lines(mapped_coord)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_triangulation.h
--- a/yt/utilities/lib/mesh_triangulation.h
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -12,7 +12,7 @@
// 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
+ {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
@@ -22,7 +22,7 @@
// Similarly, this is used to triangulate the tetrahedral cells
int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 3},
+ {0, 1, 3},
{2, 3, 1},
{0, 3, 2},
{0, 2, 1},
@@ -57,10 +57,18 @@
// 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},
+ {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}
};
+
+// This is used to select faces from a second-order tet element
+int tet10_faces[4][6] = {
+ {0, 1, 3, 4, 8, 7},
+ {2, 3, 1, 9, 8, 5},
+ {0, 3, 2, 7, 9, 6},
+ {0, 2, 1, 6, 5, 4}
+};
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_triangulation.pyx
--- a/yt/utilities/lib/mesh_triangulation.pyx
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef struct TriNode:
np.uint64_t key
@@ -35,7 +36,7 @@
if not found:
return 0
return 1
-
+
cdef np.uint64_t triangle_hash(np.int64_t tri[3]) nogil:
# http://stackoverflow.com/questions/1536393/good-hash-function-for-permutations
cdef np.uint64_t h = 1
@@ -58,13 +59,13 @@
cdef TriNode **table
cdef np.uint64_t num_items
-
+
def __cinit__(self):
self.table = <TriNode**> malloc(TABLE_SIZE * sizeof(TriNode*))
for i in range(TABLE_SIZE):
self.table[i] = NULL
self.num_items = 0
-
+
def __dealloc__(self):
cdef np.int64_t i
cdef TriNode *node
@@ -77,7 +78,7 @@
free(delete_node)
self.table[i] = NULL
free(self.table)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_exterior_tris(self):
@@ -102,7 +103,7 @@
element_map[counter] = node.elem
counter += 1
node = node.next_node
-
+
return tri_indices, element_map
cdef TriNode* _allocate_new_node(self,
@@ -118,13 +119,13 @@
new_node.next_node = NULL
self.num_items += 1
return new_node
-
+
@cython.cdivision(True)
cdef void update(self, np.int64_t tri[3], np.int64_t elem) nogil:
cdef np.uint64_t key = triangle_hash(tri)
cdef np.uint64_t index = key % TABLE_SIZE
cdef TriNode *node = self.table[index]
-
+
if node == NULL:
self.table[index] = self._allocate_new_node(tri, key, elem)
return
@@ -139,7 +140,7 @@
elif node.next_node == NULL:
node.next_node = self._allocate_new_node(tri, key, elem)
return
-
+
# walk through node list
cdef TriNode* prev = node
node = node.next_node
@@ -165,7 +166,7 @@
cdef np.int64_t VPE # num verts per element
cdef np.int64_t TPE # num tris per element
cdef int[MAX_NUM_TRI][3] tri_array
-
+
def __cinit__(self, np.int64_t[:, ::1] indices):
'''
@@ -179,7 +180,7 @@
if (self.VPE == 8 or self.VPE == 20 or self.VPE == 27):
self.TPE = HEX_NT
self.tri_array = triangulate_hex
- elif self.VPE == 4:
+ elif (self.VPE == 4 or self.VPE == 10):
self.TPE = TETRA_NT
self.tri_array = triangulate_tetra
elif self.VPE == 6:
@@ -190,7 +191,7 @@
self.num_tri = self.TPE * self.num_elem
self.num_verts = self.num_tri * 3
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def cull_interior_triangles(np.int64_t[:, ::1] indices):
@@ -213,7 +214,7 @@
s.update(tri, i)
return s.get_exterior_tris()
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_vertex_data(np.float64_t[:, ::1] coords,
@@ -230,7 +231,7 @@
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t num_verts = coords.shape[0]
cdef np.float32_t[:] vertex_data = np.zeros(num_verts, dtype="float32")
-
+
cdef np.int64_t i, j
for i in range(m.num_elem):
for j in range(m.VPE):
@@ -256,17 +257,17 @@
cdef np.int64_t num_tri = exterior_tris.shape[0]
cdef np.int64_t num_verts = 3 * num_tri
cdef np.int64_t num_coords = 3 * num_verts
-
+
cdef np.float32_t[:] vertex_data
if data.ndim == 2:
vertex_data = get_vertex_data(coords, data, indices)
else:
vertex_data = data.astype("float32")
-
+
cdef np.int32_t[:] tri_indices = np.empty(num_verts, dtype=np.int32)
cdef np.float32_t[:] tri_data = np.empty(num_verts, dtype=np.float32)
cdef np.float32_t[:] tri_coords = np.empty(num_coords, dtype=np.float32)
-
+
cdef np.int64_t vert_index, i, j, k
for i in range(num_tri):
for j in range(3):
@@ -278,7 +279,7 @@
tri_indices[vert_index] = vert_index
for k in range(3):
tri_coords[vert_index*3 + k] = coords[exterior_tris[i, j], k]
-
+
return np.array(tri_coords), np.array(tri_data), np.array(tri_indices)
@cython.boundscheck(False)
@@ -291,10 +292,10 @@
coordinates and the data values.
'''
-
+
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t[:, ::1] tri_indices = np.empty((m.num_tri, 3), dtype=np.int64)
-
+
cdef np.int64_t i, j, k
for i in range(m.num_elem):
for j in range(m.TPE):
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/primitives.pxd
--- a/yt/utilities/lib/primitives.pxd
+++ b/yt/utilities/lib/primitives.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
cimport cython.floating
import numpy as np
cimport numpy as np
@@ -66,7 +66,7 @@
cdef RayHitData compute_patch_hit(cython.floating[8][3] verts,
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil
-
+
cdef np.int64_t ray_patch_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil
@@ -78,3 +78,38 @@
cdef void patch_bbox(const void *primitives,
const np.int64_t item,
BBox* bbox) nogil
+
+cdef struct TetPatch:
+ np.float64_t[6][3] v # 6 vertices per patch
+ np.int64_t elem_id
+
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil
+
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil
+
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil
+
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil
+
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+cdef void tet_patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+cdef void tet_patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -21,15 +21,15 @@
cdef np.float64_t tmin = -INF
cdef np.float64_t tmax = INF
-
+
cdef np.int64_t i
cdef np.float64_t t1, t2
for i in range(3):
t1 = (bbox.left_edge[i] - ray.origin[i])*ray.inv_dir[i]
- t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
+ t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
tmin = fmax(tmin, fmin(t1, t2))
tmax = fmin(tmax, fmax(t1, t2))
-
+
return tmax >= fmax(tmin, 0.0)
@cython.boundscheck(False)
@@ -53,7 +53,7 @@
cdef np.float64_t det, inv_det
det = dot(e1, P)
- if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
+ if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
return False
inv_det = 1.0 / det
@@ -87,7 +87,7 @@
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
for i in range(3):
@@ -112,7 +112,7 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceFunc(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] S) nogil:
@@ -132,9 +132,9 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivU(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
- cython.floating[3] Su) nogil:
+ cython.floating[3] Su) nogil:
cdef int i
for i in range(3):
Su[i] = (-0.25*(v - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
@@ -149,10 +149,10 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivV(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil:
- cdef int i
+ cdef int i
for i in range(3):
Sv[i] = (-0.25*(u - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
(-0.25*(u + 1.0)*(u - v - 1) + 0.25*(u + 1.0)*(v - 1.0))*verts[1][i] + \
@@ -196,7 +196,7 @@
cdef cython.floating fu = dot(N1, S) + d1
cdef cython.floating fv = dot(N2, S) + d2
cdef cython.floating err = fmax(fabs(fu), fabs(fv))
-
+
# begin Newton interation
cdef cython.floating tol = 1.0e-5
cdef int iterations = 0
@@ -213,11 +213,11 @@
J21 = dot(N2, Su)
J22 = dot(N2, Sv)
det = (J11*J22 - J12*J21)
-
+
# update the u, v values
u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
-
+
patchSurfaceFunc(verts, u, v, S)
fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
@@ -227,7 +227,7 @@
# t is the distance along the ray to this hit
cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
-
+
# return hit data
cdef RayHitData hd
hd.u = u
@@ -247,7 +247,7 @@
cdef Patch patch = (<Patch*> primitives)[item]
cdef RayHitData hd = compute_patch_hit(patch.v, ray.origin, ray.direction)
-
+
# only count this is it's the closest hit
if (hd.t < ray.t_near or hd.t > ray.t_far):
return False
@@ -259,7 +259,7 @@
return True
return False
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -291,7 +291,7 @@
cdef np.int64_t i, j
cdef Patch patch = (<Patch*> primitives)[item]
-
+
for j in range(3):
bbox.left_edge[j] = patch.v[0][j]
bbox.right_edge[j] = patch.v[0][j]
@@ -300,3 +300,196 @@
for j in range(3):
bbox.left_edge[j] = fmin(bbox.left_edge[j], patch.v[i][j])
bbox.right_edge[j] = fmax(bbox.right_edge[j], patch.v[i][j])
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ S[i] = (1.0 - 3.0*u + 2.0*u*u - 3.0*v + 2.0*v*v + 4.0*u*v)*verts[0][i] + \
+ (-u + 2.0*u*u)*verts[1][i] + \
+ (-v + 2.0*v*v)*verts[2][i] + \
+ (4.0*u - 4.0*u*u - 4.0*u*v)*verts[3][i] + \
+ (4.0*u*v)*verts[4][i] + \
+ (4.0*v - 4.0*v*v - 4.0*u*v)*verts[5][i]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil:
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Su[i] = (-3.0 + 4.0*u + 4.0*v)*verts[0][i] + \
+ (-1.0 + 4.0*u)*verts[1][i] + \
+ (4.0 - 8.0*u - 4.0*v)*verts[3][i] + \
+ (4.0*v)*verts[4][i] + \
+ (-4.0*v)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Sv[i] = (-3.0 + 4.0*v + 4.0*u)*verts[0][i] + \
+ (-1.0 + 4.0*v)*verts[2][i] + \
+ (-4.0*u)*verts[3][i] + \
+ (4.0*u)*verts[4][i] + \
+ (4.0 - 8.0*v - 4.0*u)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil:
+
+ # first we compute the two planes that define the ray.
+ cdef cython.floating[3] n, N1, N2
+ cdef cython.floating A = dot(ray_direction, ray_direction)
+ for i in range(3):
+ n[i] = ray_direction[i] / A
+
+ if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):
+ N1[0] = n[1]
+ N1[1] =-n[0]
+ N1[2] = 0.0
+ else:
+ N1[0] = 0.0
+ N1[1] = n[2]
+ N1[2] =-n[1]
+ cross(N1, n, N2)
+
+ cdef cython.floating d1 = -dot(N1, ray_origin)
+ cdef cython.floating d2 = -dot(N2, ray_origin)
+
+ # the initial guess is set to zero
+ cdef cython.floating u = 0.0
+ cdef cython.floating v = 0.0
+ cdef cython.floating[3] S
+ tet_patchSurfaceFunc(verts, u, v, S)
+ cdef cython.floating fu = dot(N1, S) + d1
+ cdef cython.floating fv = dot(N2, S) + d2
+ cdef cython.floating err = fmax(fabs(fu), fabs(fv))
+
+ # begin Newton interation
+ cdef cython.floating tol = 1.0e-5
+ cdef int iterations = 0
+ cdef int max_iter = 10
+ cdef cython.floating[3] Su
+ cdef cython.floating[3] Sv
+ cdef cython.floating J11, J12, J21, J22, det
+ while ((err > tol) and (iterations < max_iter)):
+ # compute the Jacobian
+ tet_patchSurfaceDerivU(verts, u, v, Su)
+ tet_patchSurfaceDerivV(verts, u, v, Sv)
+ J11 = dot(N1, Su)
+ J12 = dot(N1, Sv)
+ J21 = dot(N2, Su)
+ J22 = dot(N2, Sv)
+ det = (J11*J22 - J12*J21)
+
+ # update the u, v values
+ u -= ( J22*fu - J12*fv) / det
+ v -= (-J21*fu + J11*fv) / det
+
+ tet_patchSurfaceFunc(verts, u, v, S)
+ fu = dot(N1, S) + d1
+ fv = dot(N2, S) + d2
+
+ err = fmax(fabs(fu), fabs(fv))
+ iterations += 1
+
+ # t is the distance along the ray to this hit
+ cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
+
+ # return hit data
+ cdef RayHitData hd
+ hd.u = u
+ hd.v = v
+ hd.t = t
+ hd.converged = (iterations < max_iter)
+ return hd
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
+
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.origin, ray.direction)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.t_near or hd.t > ray.t_far):
+ return False
+
+ if (0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1.0 and hd.converged):
+ # we have a hit, so update ray information
+ ray.t_far = hd.t
+ ray.elem_id = tet_patch.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ centroid[j] = 0.0
+
+ for i in range(6):
+ for j in range(3):
+ centroid[j] += tet_patch.v[i][j]
+
+ for j in range(3):
+ centroid[j] /= 6.0
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ bbox.left_edge[j] = tet_patch.v[0][j]
+ bbox.right_edge[j] = tet_patch.v[0][j]
+
+ for i in range(1, 6):
+ for j in range(3):
+ bbox.left_edge[j] = fmin(bbox.left_edge[j], tet_patch.v[i][j])
+ bbox.right_edge[j] = fmax(bbox.right_edge[j], tet_patch.v[i][j])
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/orientation.py
--- a/yt/utilities/orientation.py
+++ b/yt/utilities/orientation.py
@@ -23,7 +23,7 @@
def _aligned(a, b):
aligned_component = np.abs(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
return np.isclose(aligned_component, 1.0, 1.0e-13)
-
+
def _validate_unit_vectors(normal_vector, north_vector):
@@ -35,7 +35,6 @@
if not np.dot(normal_vector, normal_vector) > 0:
raise YTException("normal_vector cannot be the zero vector.")
-
if north_vector is not None and _aligned(north_vector, normal_vector):
raise YTException("normal_vector and north_vector cannot be aligned.")
@@ -85,7 +84,7 @@
t = np.cross(normal_vector, vecs).sum(axis=1)
ax = t.argmax()
east_vector = np.cross(vecs[ax, :], normal_vector).ravel()
- # self.north_vector must remain None otherwise rotations about a fixed axis will break.
+ # self.north_vector must remain None otherwise rotations about a fixed axis will break.
# The north_vector calculated here will still be included in self.unit_vectors.
north_vector = np.cross(normal_vector, east_vector).ravel()
else:
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -158,7 +158,7 @@
def position():
doc = '''
- The location of the camera.
+ The location of the camera.
Parameters
----------
@@ -179,7 +179,7 @@
'Cannot set the camera focus and position to the same value')
self._position = position
self.switch_orientation(normal_vector=self.focus - self._position,
- north_vector=None)
+ north_vector=self.north_vector)
def fdel(self):
del self._position
@@ -386,9 +386,9 @@
calculated automatically.
"""
+ if north_vector is not None:
+ self.north_vector = north_vector
self.position = position
- self.switch_orientation(normal_vector=self.focus - self.position,
- north_vector=north_vector)
def get_position(self):
"""Return the current camera position"""
@@ -483,11 +483,11 @@
>>> sc = Scene()
>>> cam = sc.add_camera()
>>> # rotate the camera by pi / 4 radians:
- >>> cam.rotate(np.pi/4.0)
+ >>> cam.rotate(np.pi/4.0)
>>> # rotate the camera about the y-axis instead of cam.north_vector:
- >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
>>> # rotate the camera about the origin instead of its own position:
- >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
"""
rotate_all = rot_vector is not None
@@ -593,7 +593,7 @@
>>> sc = Scene()
>>> cam = sc.add_camera(ds)
>>> # roll the camera by pi / 4 radians:
- >>> cam.roll(np.pi/4.0)
+ >>> cam.roll(np.pi/4.0)
>>> # roll the camera about the origin instead of its own position:
>>> cam.roll(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
@@ -627,7 +627,7 @@
>>> import yt
>>> import numpy as np
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
>>> cam = sc.camera
>>> for i in cam.iter_rotate(np.pi, 10):
@@ -688,7 +688,7 @@
def zoom(self, factor):
r"""Change the width of the FOV of the camera.
- This will appear to zoom the camera in by some `factor` toward the
+ This will appear to zoom the camera in by some `factor` toward the
focal point along the current view direction, but really it's just
changing the width of the field of view.
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -131,7 +131,7 @@
--------
The easiest way to make a VolumeSource is to use the volume_render
- function, so that the VolumeSource gets created automatically. This
+ function, so that the VolumeSource gets created automatically. This
example shows how to do this and then access the resulting source:
>>> import yt
@@ -545,7 +545,7 @@
'''
This is the name of the colormap that will be used when rendering
this MeshSource object. Should be a string, like 'arbre', or 'dusk'.
-
+
'''
def fget(self):
@@ -562,7 +562,7 @@
'''
These are the bounds that will be used with the colormap to the display
the rendered image. Should be a (vmin, vmax) tuple, like (0.0, 2.0). If
- None, the bounds will be automatically inferred from the max and min of
+ None, the bounds will be automatically inferred from the max and min of
the rendered data.
'''
@@ -603,10 +603,10 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry. Right now, high-order geometry is only
# implemented for 20-point hexes.
- if indices.shape[1] == 20:
+ if indices.shape[1] == 20 or indices.shape[1] == 10:
self.mesh = mesh_construction.QuadraticElementMesh(self.volume,
vertices,
indices,
@@ -620,12 +620,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
self.mesh = mesh_construction.LinearElementMesh(self.volume,
vertices,
@@ -651,7 +645,7 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry.
if indices.shape[1] == 27:
# hexahedral
@@ -659,12 +653,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
self.volume = BVH(vertices, indices, field_data)
@@ -794,7 +782,7 @@
cmap_name=self._cmap)/255.
alpha = image[:, :, 3]
alpha[self.sampler.aimage_used == -1] = 0.0
- image[:, :, 3] = alpha
+ image[:, :, 3] = alpha
return image
def __repr__(self):
@@ -854,7 +842,7 @@
assert(positions.ndim == 2 and positions.shape[1] == 3)
if colors is not None:
assert(colors.ndim == 2 and colors.shape[1] == 4)
- assert(colors.shape[0] == positions.shape[0])
+ assert(colors.shape[0] == positions.shape[0])
self.positions = positions
# If colors aren't individually set, make black with full opacity
if colors is None:
@@ -1057,7 +1045,7 @@
Examples
--------
- This example shows how to use BoxSource to add an outline of the
+ This example shows how to use BoxSource to add an outline of the
domain boundaries to a volume rendering.
>>> import yt
@@ -1065,12 +1053,12 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> box_source = BoxSource(ds.domain_left_edge,
... ds.domain_right_edge,
... [1.0, 1.0, 1.0, 1.0])
>>> sc.add_source(box_source)
- >>>
+ >>>
>>> im = sc.render()
"""
@@ -1078,7 +1066,7 @@
assert(left_edge.shape == (3,))
assert(right_edge.shape == (3,))
-
+
if color is None:
color = np.array([1.0, 1.0, 1.0, 1.0])
@@ -1118,7 +1106,7 @@
Examples
--------
- This example makes a volume rendering and adds outlines of all the
+ This example makes a volume rendering and adds outlines of all the
AMR grids in the simulation:
>>> import yt
@@ -1140,12 +1128,12 @@
>>> import yt
>>> from yt.visualization.volume_rendering.api import GridSource
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> dd = ds.sphere("c", (0.1, "unitary"))
>>> grid_source = GridSource(dd, alpha=1.0)
- >>>
+ >>>
>>> sc.add_source(grid_source)
>>>
>>> im = sc.render()
@@ -1223,11 +1211,11 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> coord_source = CoordinateVectorSource()
- >>>
+ >>>
>>> sc.add_source(coord_source)
- >>>
+ >>>
>>> im = sc.render()
"""
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/tests/test_camera_attributes.py
--- a/yt/visualization/volume_rendering/tests/test_camera_attributes.py
+++ b/yt/visualization/volume_rendering/tests/test_camera_attributes.py
@@ -110,3 +110,12 @@
for lens_type in valid_lens_types:
cam.set_lens(lens_type)
+
+ # See issue #1287
+ cam.focus = [0, 0, 0]
+ cam_pos = [1, 0, 0]
+ north_vector = [0, 1, 0]
+ cam.set_position(cam_pos, north_vector)
+ cam_pos = [0, 1, 0]
+ north_vector = [0, 0, 1]
+ cam.set_position(cam_pos, north_vector)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/tests/test_mesh_render.py
--- a/yt/visualization/volume_rendering/tests/test_mesh_render.py
+++ b/yt/visualization/volume_rendering/tests/test_mesh_render.py
@@ -55,7 +55,7 @@
images.append(im)
return images
-
+
@requires_module("pyembree")
def test_surface_mesh_render_pyembree():
ytcfg["yt", "ray_tracing_engine"] = "embree"
@@ -145,7 +145,7 @@
im = sc.render()
return compare(ds, im, "%s_render_answers_wedge6_%s_%s" %
(engine, field[0], field[1]))
-
+
@requires_ds(wedge6)
@requires_module("pyembree")
def test_wedge6_render_pyembree():
@@ -157,6 +157,28 @@
for field in wedge6_fields:
yield wedge6_render("yt", field)
+tet10 = "SecondOrderTets/tet10_unstructured_out.e"
+tet10_fields = [('connect1', 'uz')]
+
+def tet10_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(tet10, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_tet10_%s_%s" %
+ (engine, field[0], field[1]))
+
+ at requires_ds(tet10)
+ at requires_module("pyembree")
+def test_tet10_render_pyembree():
+ for field in tet10_fields:
+ yield tet10_render("embree", field)
+
+ at requires_ds(tet10)
+def test_tet10_render():
+ for field in tet10_fields:
+ yield tet10_render("yt", field)
+
def perspective_mesh_render(engine):
ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(hex8)
@@ -169,7 +191,7 @@
cam.resolution = (800, 800)
im = sc.render()
return compare(ds, im, "%s_perspective_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_perspective_mesh_render_pyembree():
@@ -195,7 +217,7 @@
sc.add_source(ms2)
im = sc.render()
return compare(ds, im, "%s_composite_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_composite_mesh_render_pyembree():
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