[yt-svn] commit/yt: xarthisius: Merged in atmyers/yt (pull request #2143)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed May 11 13:04:18 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/6068831a37f4/
Changeset: 6068831a37f4
Branch: yt
User: xarthisius
Date: 2016-05-11 20:04:07+00:00
Summary: Merged in atmyers/yt (pull request #2143)
High-Order Hexes for the Cython Ray Tracer
Affected #: 15 files
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 setup.py
--- a/setup.py
+++ b/setup.py
@@ -157,14 +157,21 @@
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),
+ Extension("yt.utilities.lib.primitives",
+ ["yt/utilities/lib/primitives.pyx"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/primitives.pxd",
+ "yt/utilities/lib/vec3_ops.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.bounding_volume_hierarchy",
["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
include_dirs=["yt/utilities/lib/"],
extra_compile_args=omp_args,
extra_link_args=omp_args,
libraries=std_libs,
- depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",
- "yt/utilities/lib/vec3_ops.pxd"]),
+ depends=["yt/utilities/lib/element_mappings.pxd",
+ "yt/utilities/lib/vec3_ops.pxd",
+ "yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],
include_dirs=["yt/utilities/lib/",
@@ -275,20 +282,30 @@
embree_extensions = [
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],
- depends=["yt/utilities/lib/mesh_construction.pxd"]),
+ depends=["yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/mesh_intersection.pxd",
+ "yt/utlilites/lib/mesh_samplers.pxd",
+ "yt/utlilites/lib/mesh_traversal.pxd"]),
Extension("yt.utilities.lib.mesh_traversal",
["yt/utilities/lib/mesh_traversal.pyx"],
depends=["yt/utilities/lib/mesh_traversal.pxd",
- "yt/utilities/lib/grid_traversal.pxd"]),
+ "yt/utilities/lib/grid_traversal.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.mesh_samplers",
["yt/utilities/lib/mesh_samplers.pyx"],
depends=["yt/utilities/lib/mesh_samplers.pxd",
"yt/utilities/lib/element_mappings.pxd",
- "yt/utilities/lib/mesh_construction.pxd"]),
+ "yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd",
+ "yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.mesh_intersection",
["yt/utilities/lib/mesh_intersection.pyx"],
depends=["yt/utilities/lib/mesh_intersection.pxd",
- "yt/utilities/lib/mesh_construction.pxd"]),
+ "yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd",
+ "yt/utilities/lib/mesh_samplers.pxd",
+ "yt/utilities/lib/primitives.pxd",
+ "yt/utilities/lib/vec3_ops.pxd"]),
]
embree_prefix = os.path.abspath(read_embree_location())
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -47,7 +47,7 @@
local_tipsy_000:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_000:
+ local_varia_001:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/config.py
--- a/yt/config.py
+++ b/yt/config.py
@@ -65,6 +65,7 @@
chunk_size = '1000',
xray_data_dir = '/does/not/exist',
default_colormap = 'arbre',
+ ray_tracing_engine = 'embree',
)
# Here is the upgrade. We're actually going to parse the file in its entirety
# here. Then, if it has any of the Forbidden Sections, it will be rewritten
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/testing.py
--- a/yt/testing.py
+++ b/yt/testing.py
@@ -306,18 +306,16 @@
from yt.frontends.stream.sample_data.hexahedral_mesh import \
_connectivity, _coordinates
- _connectivity -= 1 # this mesh has an offset of 1
-
# the distance from the origin
node_data = {}
dist = np.sum(_coordinates**2, 1)
- node_data[('connect1', 'test')] = dist[_connectivity]
+ node_data[('connect1', 'test')] = dist[_connectivity-1]
# each element gets a random number
elem_data = {}
elem_data[('connect1', 'elem')] = np.random.rand(_connectivity.shape[0])
- ds = load_unstructured_mesh(_connectivity,
+ ds = load_unstructured_mesh(_connectivity-1,
_coordinates,
node_data=node_data,
elem_data=elem_data)
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -2,22 +2,22 @@
import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler
+from yt.utilities.lib.primitives cimport BBox, Ray
-# ray data structure
-cdef struct Ray:
- np.float64_t origin[3]
- np.float64_t direction[3]
- np.float64_t inv_dir[3]
- np.float64_t data_val
- np.float64_t t_near
- np.float64_t t_far
- np.int64_t elem_id
- np.int64_t near_boundary
+cdef extern from "mesh_construction.h":
+ enum:
+ MAX_NUM_TRI
-# axis-aligned bounding box
-cdef struct BBox:
- np.float64_t left_edge[3]
- np.float64_t right_edge[3]
+ int HEX_NV
+ int HEX_NT
+ int TETRA_NV
+ int TETRA_NT
+ int WEDGE_NV
+ int WEDGE_NT
+ int triangulate_hex[MAX_NUM_TRI][3]
+ int triangulate_tetra[MAX_NUM_TRI][3]
+ int triangulate_wedge[MAX_NUM_TRI][3]
+ int hex20_faces[6][8]
# node for the bounding volume hierarchy
cdef struct BVHNode:
@@ -26,29 +26,49 @@
BVHNode* left
BVHNode* right
BBox bbox
-
-# triangle data structure
-cdef struct Triangle:
- np.float64_t p0[3]
- np.float64_t p1[3]
- np.float64_t p2[3]
- np.int64_t elem_id
- np.float64_t centroid[3]
- BBox bbox
+
+# pointer to function that computes primitive intersection
+ctypedef np.int64_t (*intersect_func_type)(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+# pointer to function that computes primitive centroids
+ctypedef void (*centroid_func_type)(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+# pointer to function that computes primitive bounding boxes
+ctypedef void (*bbox_func_type)(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
+
cdef class BVH:
cdef BVHNode* root
- cdef Triangle* triangles
+ cdef void* primitives
+ cdef np.int64_t* prim_ids
+ cdef np.float64_t** centroids
+ cdef BBox* bboxes
cdef np.float64_t* vertices
cdef np.float64_t* field_data
- cdef np.int64_t num_tri_per_elem
- cdef np.int64_t num_tri
+ cdef np.int64_t num_prim_per_elem
+ cdef np.int64_t num_prim
cdef np.int64_t num_elem
cdef np.int64_t num_verts_per_elem
cdef np.int64_t num_field_per_elem
+ cdef int[MAX_NUM_TRI][3] tri_array
cdef ElementSampler sampler
+ cdef centroid_func_type get_centroid
+ cdef bbox_func_type get_bbox
+ cdef intersect_func_type get_intersect
cdef np.int64_t _partition(self, np.int64_t begin, np.int64_t end,
np.int64_t ax, np.float64_t split) nogil
+ cdef void _set_up_triangles(self,
+ np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil
+ cdef void _set_up_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 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -4,106 +4,41 @@
from libc.math cimport fabs
from libc.stdlib cimport malloc, free
from cython.parallel import parallel, prange
-from vec3_ops cimport dot, subtract, cross
+
+from yt.utilities.lib.primitives cimport \
+ BBox, \
+ Ray, \
+ ray_bbox_intersect, \
+ Triangle, \
+ ray_triangle_intersect, \
+ triangle_centroid, \
+ triangle_bbox, \
+ Patch, \
+ ray_patch_intersect, \
+ patch_centroid, \
+ patch_bbox
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
P1Sampler3D, \
- W1Sampler3D
+ W1Sampler3D, \
+ S2Sampler3D
+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 extern from "platform_dep.h" nogil:
double fmax(double x, double y)
double fmin(double x, double y)
-cdef extern from "mesh_construction.h":
- enum:
- MAX_NUM_TRI
-
- int HEX_NV
- int HEX_NT
- int TETRA_NV
- int TETRA_NT
- int WEDGE_NV
- int WEDGE_NT
- int triangulate_hex[MAX_NUM_TRI][3]
- int triangulate_tetra[MAX_NUM_TRI][3]
- int triangulate_wedge[MAX_NUM_TRI][3]
-
# define some constants
-cdef np.float64_t DETERMINANT_EPS = 1.0e-10
cdef np.float64_t INF = np.inf
cdef np.int64_t LEAF_SIZE = 16
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef np.int64_t ray_triangle_intersect(Ray* ray, const Triangle* tri) nogil:
-# https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
-
- # edge vectors
- cdef np.float64_t e1[3]
- cdef np.float64_t e2[3]
- subtract(tri.p1, tri.p0, e1)
- subtract(tri.p2, tri.p0, e2)
-
- cdef np.float64_t P[3]
- cross(ray.direction, e2, P)
-
- cdef np.float64_t det, inv_det
- det = dot(e1, P)
- if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
- return False
- inv_det = 1.0 / det
-
- cdef np.float64_t T[3]
- subtract(ray.origin, tri.p0, T)
-
- cdef np.float64_t u = dot(T, P) * inv_det
- if(u < 0.0 or u > 1.0):
- return False
-
- cdef np.float64_t Q[3]
- cross(T, e1, Q)
-
- cdef np.float64_t v = dot(ray.direction, Q) * inv_det
- if(v < 0.0 or u + v > 1.0):
- return False
-
- cdef np.float64_t t = dot(e2, Q) * inv_det
-
- if(t > DETERMINANT_EPS and t < ray.t_far):
- ray.t_far = t
- ray.elem_id = tri.elem_id
- return True
-
- return False
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil:
-# https://tavianator.com/fast-branchless-raybounding-box-intersections/
-
- 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]
- tmin = fmax(tmin, fmin(t1, t2))
- tmax = fmin(tmax, fmax(t1, t2))
-
- return tmax >= fmax(tmin, 0.0)
-
-
cdef class BVH:
'''
@@ -130,29 +65,40 @@
self.num_field_per_elem = field_data.shape[1]
# We need to figure out what kind of elements we've been handed.
- cdef int[MAX_NUM_TRI][3] tri_array
if self.num_verts_per_elem == 8:
- self.num_tri_per_elem = HEX_NT
- tri_array = triangulate_hex
+ self.num_prim_per_elem = HEX_NT
+ self.tri_array = triangulate_hex
self.sampler = Q1Sampler
elif self.num_verts_per_elem == 6:
- self.num_tri_per_elem = WEDGE_NT
- tri_array = triangulate_wedge
+ self.num_prim_per_elem = WEDGE_NT
+ self.tri_array = triangulate_wedge
self.sampler = W1Sampler
elif self.num_verts_per_elem == 4:
- self.num_tri_per_elem = TETRA_NT
- tri_array = triangulate_tetra
+ self.num_prim_per_elem = TETRA_NT
+ self.tri_array = triangulate_tetra
self.sampler = P1Sampler
- self.num_tri = self.num_tri_per_elem*self.num_elem
+ elif self.num_verts_per_elem == 20:
+ self.num_prim_per_elem = 6
+ self.sampler = S2Sampler
+ else:
+ raise NotImplementedError("Could not determine element type for "
+ "nverts = %d. " % self.num_verts_per_elem)
+ self.num_prim = self.num_prim_per_elem*self.num_elem
# allocate storage
cdef np.int64_t v_size = self.num_verts_per_elem * self.num_elem * 3
self.vertices = <np.float64_t*> malloc(v_size * sizeof(np.float64_t))
cdef np.int64_t f_size = self.num_field_per_elem * self.num_elem
self.field_data = <np.float64_t*> malloc(f_size * sizeof(np.float64_t))
+ self.prim_ids = <np.int64_t*> malloc(self.num_prim * sizeof(np.int64_t))
+ self.centroids = <np.float64_t**> malloc(self.num_prim * sizeof(np.float64_t*))
+ cdef np.int64_t i
+ for i in range(self.num_prim):
+ self.centroids[i] = <np.float64_t*> malloc(3*sizeof(np.float64_t))
+ self.bboxes = <BBox*> malloc(self.num_prim * sizeof(BBox))
# create data buffers
- cdef np.int64_t i, j, k
+ cdef np.int64_t j, k
cdef np.int64_t field_offset, vertex_offset
for i in range(self.num_elem):
for j in range(self.num_verts_per_elem):
@@ -163,29 +109,78 @@
for j in range(self.num_field_per_elem):
self.field_data[field_offset + j] = field_data[i][j]
- # fill our array of triangles
+ # set up primitives
+ if self.num_verts_per_elem == 20:
+ self.primitives = malloc(self.num_prim * sizeof(Patch))
+ self.get_centroid = patch_centroid
+ self.get_bbox = patch_bbox
+ self.get_intersect = ray_patch_intersect
+ self._set_up_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)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void _set_up_patches(self, np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil:
+ cdef Patch* 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
+ patch = &( <Patch*> self.primitives)[prim_index]
+ self.prim_ids[prim_index] = prim_index
+ patch.elem_id = i
+ for k in range(8): # for each vertex
+ ind = hex20_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ 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
cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri
- self.triangles = <Triangle*> malloc(self.num_tri * sizeof(Triangle))
+ cdef np.int64_t i, j, k
for i in range(self.num_elem):
- offset = self.num_tri_per_elem*i
- for j in range(self.num_tri_per_elem):
+ offset = self.num_prim_per_elem*i
+ for j in range(self.num_prim_per_elem):
tri_index = offset + j
- tri = &(self.triangles[tri_index])
+ self.prim_ids[tri_index] = tri_index
+ tri = &(<Triangle*> self.primitives)[tri_index]
tri.elem_id = i
- v0 = indices[i][tri_array[j][0]]
- v1 = indices[i][tri_array[j][1]]
- v2 = indices[i][tri_array[j][2]]
+ v0 = indices[i][self.tri_array[j][0]]
+ v1 = indices[i][self.tri_array[j][1]]
+ v2 = indices[i][self.tri_array[j][2]]
for k in range(3):
tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]
- tri.centroid[k] = (tri.p0[k] + tri.p1[k] + tri.p2[k]) / 3.0
- tri.bbox.left_edge[k] = fmin(fmin(tri.p0[k], tri.p1[k]), tri.p2[k])
- tri.bbox.right_edge[k] = fmax(fmax(tri.p0[k], tri.p1[k]), tri.p2[k])
-
- self.root = self._recursive_build(0, self.num_tri)
+ self.get_centroid(self.primitives,
+ tri_index,
+ self.centroids[tri_index])
+ self.get_bbox(self.primitives,
+ tri_index,
+ &(self.bboxes[tri_index]))
cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:
@@ -195,7 +190,12 @@
def __dealloc__(self):
self._recursive_free(self.root)
- free(self.triangles)
+ free(self.primitives)
+ free(self.prim_ids)
+ for i in range(self.num_prim):
+ free(self.centroids[i])
+ free(self.centroids)
+ free(self.bboxes)
free(self.field_data)
free(self.vertices)
@@ -204,17 +204,21 @@
@cython.cdivision(True)
cdef np.int64_t _partition(self, np.int64_t begin, np.int64_t end,
np.int64_t ax, np.float64_t split) nogil:
- # this re-orders the triangle array so that all of the triangles
- # to the left of mid have centroids less than or equal to "split"
- # along the direction "ax". All the triangles to the right of mid
+ # this re-orders the primitive array so that all of the primitives
+ # to the left of mid have centroids less than or equal to "split"
+ # along the direction "ax". All the primitives to the right of mid
# will have centroids *greater* than "split" along "ax".
cdef np.int64_t mid = begin
while (begin != end):
- if self.triangles[mid].centroid[ax] > split:
+ if self.centroids[mid][ax] > split:
mid += 1
- elif self.triangles[begin].centroid[ax] > split:
- self.triangles[mid], self.triangles[begin] = \
- self.triangles[begin], self.triangles[mid]
+ elif self.centroids[begin][ax] > split:
+ self.prim_ids[mid], self.prim_ids[begin] = \
+ self.prim_ids[begin], self.prim_ids[mid]
+ self.centroids[mid], self.centroids[begin] = \
+ self.centroids[begin], self.centroids[mid]
+ self.bboxes[mid], self.bboxes[begin] = \
+ self.bboxes[begin], self.bboxes[mid]
mid += 1
begin += 1
return mid
@@ -225,13 +229,13 @@
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.triangles[begin].bbox
+ cdef BBox box = self.bboxes[begin]
for i in range(begin+1, end):
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],
- self.triangles[i].bbox.left_edge[j])
+ self.bboxes[i].left_edge[j])
box.right_edge[j] = fmax(box.right_edge[j],
- self.triangles[i].bbox.right_edge[j])
+ self.bboxes[i].right_edge[j])
node.bbox = box
@cython.boundscheck(False)
@@ -249,7 +253,7 @@
position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
cdef np.float64_t* vertex_ptr
- cdef np.float64_t* field_ptr
+ cdef np.float64_t* field_ptr
vertex_ptr = self.vertices + ray.elem_id*self.num_verts_per_elem*3
field_ptr = self.field_data + ray.elem_id*self.num_field_per_elem
@@ -273,11 +277,9 @@
# check for leaf
cdef np.int64_t i, hit
- cdef Triangle* tri
if (node.end - node.begin) <= LEAF_SIZE:
for i in range(node.begin, node.end):
- tri = &(self.triangles[i])
- hit = ray_triangle_intersect(ray, tri)
+ hit = self.get_intersect(self.primitives, self.prim_ids[i], ray)
return
# if not leaf, intersect with left and right children
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -17,25 +17,27 @@
import numpy as np
cimport cython
+from libc.stdlib cimport malloc, free
+from libc.math cimport fmax, sqrt
+cimport numpy as np
+
cimport pyembree.rtcore as rtc
-from mesh_traversal cimport YTEmbreeScene
cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu
+from pyembree.rtcore cimport \
+ Vertex, \
+ Triangle, \
+ Vec3f
+
+from mesh_traversal cimport YTEmbreeScene
from mesh_samplers cimport \
sample_hex, \
sample_tetra, \
sample_wedge
-from pyembree.rtcore cimport \
- Vertex, \
- Triangle, \
- Vec3f
from mesh_intersection cimport \
patchIntersectFunc, \
patchBoundsFunc
-from libc.stdlib cimport malloc, free
-from libc.math cimport fmax, sqrt
-cimport numpy as np
from yt.utilities.exceptions import YTElementTypeNotRecognized
cdef extern from "mesh_construction.h":
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/mesh_intersection.pxd
--- a/yt/utilities/lib/mesh_intersection.pxd
+++ b/yt/utilities/lib/mesh_intersection.pxd
@@ -4,15 +4,10 @@
from yt.utilities.lib.mesh_construction cimport Patch
cimport cython
-cdef void patchIntersectFunc(Patch* patches,
- rtcr.RTCRay& ray,
+cdef void patchIntersectFunc(Patch* patches,
+ rtcr.RTCRay& ray,
size_t item) nogil
-cdef void patchBoundsFunc(Patch* patches,
- size_t item,
+cdef void patchBoundsFunc(Patch* patches,
+ size_t item,
rtcg.RTCBounds* bounds_o) nogil
-
-cdef void patchSurfaceFunc(const Patch& patch,
- const float u,
- const float v,
- float[3] S) nogil
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/mesh_intersection.pyx
--- a/yt/utilities/lib/mesh_intersection.pyx
+++ b/yt/utilities/lib/mesh_intersection.pyx
@@ -21,66 +21,19 @@
cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
from yt.utilities.lib.mesh_samplers cimport sample_hex20
+from yt.utilities.lib.bounding_volume_hierarchy cimport BBox
+from yt.utilities.lib.primitives cimport \
+ patchSurfaceFunc, \
+ patchSurfaceDerivU, \
+ patchSurfaceDerivV, \
+ RayHitData, \
+ compute_patch_hit
from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void patchSurfaceFunc(const Patch& patch,
- const float u,
- const float v,
- float[3] S) nogil:
-
- cdef int i
- for i in range(3):
- S[i] = 0.25*(1.0 - u)*(1.0 - v)*(-u - v - 1)*patch.v[0][i] + \
- 0.25*(1.0 + u)*(1.0 - v)*( u - v - 1)*patch.v[1][i] + \
- 0.25*(1.0 + u)*(1.0 + v)*( u + v - 1)*patch.v[2][i] + \
- 0.25*(1.0 - u)*(1.0 + v)*(-u + v - 1)*patch.v[3][i] + \
- 0.5*(1 - u)*(1 - v*v)*patch.v[4][i] + \
- 0.5*(1 - u*u)*(1 - v)*patch.v[5][i] + \
- 0.5*(1 + u)*(1 - v*v)*patch.v[6][i] + \
- 0.5*(1 - u*u)*(1 + v)*patch.v[7][i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void patchSurfaceDerivU(const Patch& patch,
- const float u,
- const float v,
- float[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))*patch.v[0][i] + \
- (-0.25*(v - 1.0)*(u - v - 1) - 0.25*(u + 1.0)*(v - 1.0))*patch.v[1][i] + \
- ( 0.25*(v + 1.0)*(u + v - 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \
- ( 0.25*(v + 1.0)*(u - v + 1) + 0.25*(u - 1.0)*(v + 1.0))*patch.v[3][i] + \
- 0.5*(v*v - 1.0)*patch.v[4][i] + u*(v - 1.0)*patch.v[5][i] - \
- 0.5*(v*v - 1.0)*patch.v[6][i] - u*(v + 1.0)*patch.v[7][i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void patchSurfaceDerivV(const Patch& patch,
- const float u,
- const float v,
- float[3] Sv) nogil:
- 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))*patch.v[0][i] + \
- (-0.25*(u + 1.0)*(u - v - 1) + 0.25*(u + 1.0)*(v - 1.0))*patch.v[1][i] + \
- ( 0.25*(u + 1.0)*(u + v - 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \
- ( 0.25*(u - 1.0)*(u - v + 1) - 0.25*(u - 1.0)*(v + 1.0))*patch.v[3][i] + \
- 0.5*(u*u - 1.0)*patch.v[5][i] + v*(u - 1.0)*patch.v[4][i] - \
- 0.5*(u*u - 1.0)*patch.v[7][i] - v*(u + 1.0)*patch.v[6][i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:
@@ -121,77 +74,20 @@
cdef Patch patch = patches[item]
- # first we compute the two planes that define the ray.
- cdef float[3] n, N1, N2
- cdef float A = dot(ray.dir, ray.dir)
- for i in range(3):
- n[i] = ray.dir[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 float d1 = -dot(N1, ray.org)
- cdef float d2 = -dot(N2, ray.org)
-
- # the initial guess is set to zero
- cdef float u = 0.0
- cdef float v = 0.0
- cdef float[3] S
- patchSurfaceFunc(patch, u, v, S)
- cdef float fu = dot(N1, S) + d1
- cdef float fv = dot(N2, S) + d2
- cdef float err = fmax(fabs(fu), fabs(fv))
-
- # begin Newton interation
- cdef float tol = 1.0e-5
- cdef int iterations = 0
- cdef int max_iter = 10
- cdef float[3] Su
- cdef float[3] Sv
- cdef float J11, J12, J21, J22, det
- while ((err > tol) and (iterations < max_iter)):
- # compute the Jacobian
- patchSurfaceDerivU(patch, u, v, Su)
- patchSurfaceDerivV(patch, 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
-
- patchSurfaceFunc(patch, 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 float t = distance(S, ray.org)
+ cdef RayHitData hd = compute_patch_hit(patch.v, ray.org, ray.dir)
# only count this is it's the closest hit
- if (t < ray.tnear or t > ray.Ng[0]):
+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):
return
- if (fabs(u) <= 1.0 and fabs(v) <= 1.0 and iterations < max_iter):
+ if (fabs(hd.u) <= 1.0 and fabs(hd.v) <= 1.0 and hd.converged):
# we have a hit, so update ray information
- ray.u = u
- ray.v = v
+ ray.u = hd.u
+ ray.v = hd.v
ray.geomID = patch.geomID
ray.primID = item
- ray.Ng[0] = t
+ ray.Ng[0] = hd.t
# sample the solution at the calculated point
sample_hex20(patches, ray)
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/mesh_samplers.pyx
--- a/yt/utilities/lib/mesh_samplers.pyx
+++ b/yt/utilities/lib/mesh_samplers.pyx
@@ -19,7 +19,7 @@
from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
Patch
-from yt.utilities.lib.mesh_intersection cimport patchSurfaceFunc
+from yt.utilities.lib.primitives cimport patchSurfaceFunc
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \
@@ -192,7 +192,7 @@
elem_id = ray_id / 6
# fills "position" with the physical position of the hit
- patchSurfaceFunc(data[ray_id], ray.u, ray.v, pos)
+ patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
for i in range(3):
position[i] = <double> pos[i]
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/primitives.pxd
--- /dev/null
+++ b/yt/utilities/lib/primitives.pxd
@@ -0,0 +1,80 @@
+cimport cython
+cimport cython.floating
+import numpy as np
+cimport numpy as np
+from vec3_ops cimport dot, subtract, cross
+
+cdef struct Ray:
+ np.float64_t origin[3]
+ np.float64_t direction[3]
+ np.float64_t inv_dir[3]
+ np.float64_t data_val
+ np.float64_t t_near
+ np.float64_t t_far
+ np.int64_t elem_id
+ np.int64_t near_boundary
+
+cdef struct BBox:
+ np.float64_t left_edge[3]
+ np.float64_t right_edge[3]
+
+cdef struct RayHitData:
+ np.float64_t u
+ np.float64_t v
+ np.float64_t t
+ np.int64_t converged
+
+cdef struct Triangle:
+ np.float64_t p0[3]
+ np.float64_t p1[3]
+ np.float64_t p2[3]
+ np.int64_t elem_id
+
+cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil
+
+cdef inline np.int64_t ray_triangle_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+cdef inline void triangle_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+cdef inline void triangle_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
+
+cdef struct Patch:
+ np.float64_t[8][3] v # 8 vertices per patch
+ np.int64_t elem_id
+
+cdef void patchSurfaceFunc(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil
+
+cdef void patchSurfaceDerivU(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil
+
+cdef void patchSurfaceDerivV(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil
+
+cdef RayHitData compute_patch_hit(cython.floating[8][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil
+
+cdef inline np.int64_t ray_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+cdef inline void patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+cdef inline void patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/utilities/lib/primitives.pyx
--- /dev/null
+++ b/yt/utilities/lib/primitives.pyx
@@ -0,0 +1,302 @@
+cimport cython
+import numpy as np
+cimport numpy as np
+cimport cython.floating
+from libc.math cimport fabs
+
+from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance, L2_norm
+
+cdef np.float64_t DETERMINANT_EPS = 1.0e-10
+cdef np.float64_t INF = np.inf
+
+cdef extern from "platform_dep.h" nogil:
+ double fmax(double x, double y)
+ double fmin(double x, double y)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil:
+# https://tavianator.com/fast-branchless-raybounding-box-intersections/
+
+ 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]
+ tmin = fmax(tmin, fmin(t1, t2))
+ tmax = fmin(tmax, fmax(t1, t2))
+
+ return tmax >= fmax(tmin, 0.0)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline np.int64_t ray_triangle_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
+# https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+
+ cdef Triangle tri = (<Triangle*> primitives)[item]
+
+ # edge vectors
+ cdef np.float64_t e1[3]
+ cdef np.float64_t e2[3]
+ subtract(tri.p1, tri.p0, e1)
+ subtract(tri.p2, tri.p0, e2)
+
+ cdef np.float64_t P[3]
+ cross(ray.direction, e2, P)
+
+ cdef np.float64_t det, inv_det
+ det = dot(e1, P)
+ if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
+ return False
+ inv_det = 1.0 / det
+
+ cdef np.float64_t T[3]
+ subtract(ray.origin, tri.p0, T)
+
+ cdef np.float64_t u = dot(T, P) * inv_det
+ if(u < 0.0 or u > 1.0):
+ return False
+
+ cdef np.float64_t Q[3]
+ cross(T, e1, Q)
+
+ cdef np.float64_t v = dot(ray.direction, Q) * inv_det
+ if(v < 0.0 or u + v > 1.0):
+ return False
+
+ cdef np.float64_t t = dot(e2, Q) * inv_det
+
+ if(t > DETERMINANT_EPS and t < ray.t_far):
+ ray.t_far = t
+ ray.elem_id = tri.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline 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):
+ centroid[i] = (tri.p0[i] + tri.p1[i] + tri.p2[i]) / 3.0
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void triangle_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
+
+ cdef Triangle tri = (<Triangle*> primitives)[item]
+ cdef np.int64_t i
+ for i in range(3):
+ bbox.left_edge[i] = fmin(fmin(tri.p0[i], tri.p1[i]), tri.p2[i])
+ bbox.right_edge[i] = fmax(fmax(tri.p0[i], tri.p1[i]), tri.p2[i])
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void patchSurfaceFunc(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil:
+
+ cdef int i
+ for i in range(3):
+ S[i] = 0.25*(1.0 - u)*(1.0 - v)*(-u - v - 1)*verts[0][i] + \
+ 0.25*(1.0 + u)*(1.0 - v)*( u - v - 1)*verts[1][i] + \
+ 0.25*(1.0 + u)*(1.0 + v)*( u + v - 1)*verts[2][i] + \
+ 0.25*(1.0 - u)*(1.0 + v)*(-u + v - 1)*verts[3][i] + \
+ 0.5*(1 - u)*(1 - v*v)*verts[4][i] + \
+ 0.5*(1 - u*u)*(1 - v)*verts[5][i] + \
+ 0.5*(1 + u)*(1 - v*v)*verts[6][i] + \
+ 0.5*(1 - u*u)*(1 + v)*verts[7][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void patchSurfaceDerivU(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ 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] + \
+ (-0.25*(v - 1.0)*(u - v - 1) - 0.25*(u + 1.0)*(v - 1.0))*verts[1][i] + \
+ ( 0.25*(v + 1.0)*(u + v - 1) + 0.25*(u + 1.0)*(v + 1.0))*verts[2][i] + \
+ ( 0.25*(v + 1.0)*(u - v + 1) + 0.25*(u - 1.0)*(v + 1.0))*verts[3][i] + \
+ 0.5*(v*v - 1.0)*verts[4][i] + u*(v - 1.0)*verts[5][i] - \
+ 0.5*(v*v - 1.0)*verts[6][i] - u*(v + 1.0)*verts[7][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void patchSurfaceDerivV(const cython.floating[8][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil:
+ 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] + \
+ ( 0.25*(u + 1.0)*(u + v - 1) + 0.25*(u + 1.0)*(v + 1.0))*verts[2][i] + \
+ ( 0.25*(u - 1.0)*(u - v + 1) - 0.25*(u - 1.0)*(v + 1.0))*verts[3][i] + \
+ 0.5*(u*u - 1.0)*verts[5][i] + v*(u - 1.0)*verts[4][i] - \
+ 0.5*(u*u - 1.0)*verts[7][i] - v*(u + 1.0)*verts[6][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef RayHitData compute_patch_hit(cython.floating[8][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
+ 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
+ patchSurfaceDerivU(verts, u, v, Su)
+ 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
+
+ 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 inline np.int64_t ray_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
+
+ 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
+
+ if (fabs(hd.u) <= 1.0 and fabs(hd.v) <= 1.0 and hd.converged):
+ # we have a hit, so update ray information
+ ray.t_far = hd.t
+ ray.elem_id = patch.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
+
+ cdef np.int64_t i, j
+ cdef Patch patch = (<Patch*> primitives)[item]
+
+ for j in range(3):
+ centroid[j] = 0.0
+
+ for i in range(8):
+ for j in range(3):
+ centroid[j] += patch.v[i][j]
+
+ for j in range(3):
+ centroid[j] /= 8.0
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
+
+ cdef np.int64_t i, j
+ cdef Patch patch = (<Patch*> primitives)[item]
+
+ for j in range(3):
+ bbox.left_edge[j] = patch.v[0][j]
+ bbox.right_edge[j] = patch.v[0][j]
+
+ for i in range(1, 8):
+ 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])
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -36,10 +36,12 @@
from yt.utilities.lib import mesh_traversal
except ImportError:
mesh_traversal = NotAModule("pyembree")
+ ytcfg["yt", "ray_tracing_engine"] = "yt"
try:
from yt.utilities.lib import mesh_construction
except ImportError:
mesh_construction = NotAModule("pyembree")
+ ytcfg["yt", "ray_tracing_engine"] = "yt"
class RenderSource(ParallelAnalysisInterface):
@@ -356,7 +358,7 @@
self.field = field
self.volume = None
self.current_image = None
- self.engine = 'embree'
+ self.engine = ytcfg.get("yt", "ray_tracing_engine")
# default color map
self._cmap = ytcfg.get("yt", "default_colormap")
@@ -374,11 +376,11 @@
if self.engine == 'embree':
self.volume = mesh_traversal.YTEmbreeScene()
self.build_volume_embree()
- elif self.engine == 'bvh':
+ elif self.engine == 'yt':
self.build_volume_bvh()
else:
raise NotImplementedError("Invalid ray-tracing engine selected. "
- "Choices are 'embree' and 'bvh'.")
+ "Choices are 'embree' and 'yt'.")
def cmap():
'''
@@ -491,15 +493,8 @@
field_data = np.expand_dims(field_data, 1)
# Here, we decide whether to render based on high-order or
- # low-order geometry. Right now, high-order geometry is only
- # supported by Embree.
- if indices.shape[1] == 20:
- # hexahedral
- mylog.warning("20-node hexes not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:8]
- indices = indices[:, 0:8]
- elif indices.shape[1] == 27:
+ # low-order geometry.
+ if indices.shape[1] == 27:
# hexahedral
mylog.warning("27-node hexes not yet supported, " +
"dropping to 1st order.")
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 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
@@ -23,30 +23,8 @@
MeshSource, \
Scene, \
create_scene
-
-
- at requires_module("pyembree")
-def test_surface_mesh_render():
-
- images = []
-
- ds = fake_tetrahedral_ds()
- sc = Scene()
- for field in ds.field_list:
- sc.add_source(MeshSource(ds, field))
- sc.add_camera()
- im = sc.render()
- images.append(im)
-
- ds = fake_hexahedral_ds()
- for field in ds.field_list:
- sc.add_source(MeshSource(ds, field))
- sc.add_camera()
- im = sc.render()
- images.append(im)
-
- return images
-
+from yt.config import \
+ ytcfg
def compare(ds, im, test_prefix, decimals=12):
def mesh_render_image_func(filename_prefix):
@@ -57,64 +35,132 @@
test.prefix = test_prefix
return test
+def surface_mesh_render():
+ images = []
+
+ ds = fake_tetrahedral_ds()
+ for field in ds.field_list:
+ sc = Scene()
+ sc.add_source(MeshSource(ds, field))
+ sc.add_camera()
+ im = sc.render()
+ images.append(im)
+
+ ds = fake_hexahedral_ds()
+ for field in ds.field_list:
+ sc = Scene()
+ sc.add_source(MeshSource(ds, field))
+ sc.add_camera()
+ im = sc.render()
+ images.append(im)
+
+ return images
+
+ at requires_module("pyembree")
+def test_surface_mesh_render_pyembree():
+ ytcfg["yt", "ray_tracing_engine"] = "embree"
+ surface_mesh_render()
+
+def test_surface_mesh_render():
+ ytcfg["yt", "ray_tracing_engine"] = "yt"
+ surface_mesh_render()
+
+
hex8 = "MOOSE_sample_data/out.e-s010"
hex8_fields = [('connect1', 'diffused'), ('connect2', 'convected')]
+def hex8_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(hex8, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_hex8_%s_%s" %
+ (engine, field[0], field[1]))
+
@requires_ds(hex8)
@requires_module("pyembree")
+def test_hex8_render_pyembree():
+ for field in hex8_fields:
+ yield hex8_render("embree", field)
+
+ at requires_ds(hex8)
def test_hex8_render():
for field in hex8_fields:
- ds = data_dir_load(hex8, kwargs={'step':-1})
- sc = create_scene(ds, field)
- im = sc.render()
- yield compare(ds, im, "render_answers_hex8_%s_%s" % field)
+ yield hex8_render("yt", field)
tet4 = "MOOSE_sample_data/high_order_elems_tet4_refine_out.e"
tet4_fields = [("connect1", "u")]
+def tet4_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(tet4, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_tet4_%s_%s" %
+ (engine, field[0], field[1]))
+
@requires_ds(tet4)
@requires_module("pyembree")
+def test_tet4_render_pyembree():
+ for field in tet4_fields:
+ yield tet4_render("embree", field)
+
+ at requires_ds(tet4)
def test_tet4_render():
for field in tet4_fields:
- ds = data_dir_load(tet4, kwargs={'step':-1})
- sc = create_scene(ds, field)
- im = sc.render()
- yield compare(ds, im, "render_answers_tet4_%s_%s" % field)
+ yield tet4_render("yt", field)
hex20 = "MOOSE_sample_data/mps_out.e"
hex20_fields = [('connect2', 'temp')]
+def hex20_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(hex20, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_hex20_%s_%s" %
+ (engine, field[0], field[1]))
+
@requires_ds(hex20)
@requires_module("pyembree")
+def test_hex20_render_pyembree():
+ for field in hex20_fields:
+ yield hex20_render("embree", field)
+
+ at requires_ds(hex20)
def test_hex20_render():
for field in hex20_fields:
- ds = data_dir_load(hex20, kwargs={'step':-1})
- sc = create_scene(ds, field)
- im = sc.render()
- yield compare(ds, im, "render_answers_hex20_%s_%s" % field)
+ yield hex20_render("yt", field)
wedge6 = "MOOSE_sample_data/wedge_out.e"
wedge6_fields = [('connect1', 'diffused')]
+def wedge6_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(wedge6, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ 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():
+ for field in wedge6_fields:
+ yield wedge6_render("embree", field)
+
+ at requires_ds(wedge6)
def test_wedge6_render():
for field in wedge6_fields:
- ds = data_dir_load(wedge6, kwargs={'step':-1})
- sc = create_scene(ds, field)
- im = sc.render()
- yield compare(ds, im, "render_answers_wedge6_%s_%s" % field)
+ yield wedge6_render("yt", field)
-
- at requires_ds(hex8)
- at requires_module("pyembree")
-def test_perspective_mesh_render():
+def perspective_mesh_render(engine):
+ ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(hex8)
sc = create_scene(ds, ("connect2", "diffused"))
-
cam = sc.add_camera(ds, lens_type='perspective')
cam.focus = ds.arr([0.0, 0.0, 0.0], 'code_length')
cam_pos = ds.arr([-4.5, 4.5, -4.5], 'code_length')
@@ -122,12 +168,19 @@
cam.set_position(cam_pos, north_vector)
cam.resolution = (800, 800)
im = sc.render()
- yield compare(ds, im, "perspective_mesh_render")
-
+ return compare(ds, im, "%s_perspective_mesh_render" % engine)
+
+ at requires_ds(hex8)
+ at requires_module("pyembree")
+def test_perspective_mesh_render_pyembree():
+ yield perspective_mesh_render("embree")
@requires_ds(hex8)
- at requires_module("pyembree")
-def test_composite_mesh_render():
+def test_perspective_mesh_render():
+ yield perspective_mesh_render("yt")
+
+def composite_mesh_render(engine):
+ ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(hex8)
sc = Scene()
cam = sc.add_camera(ds)
@@ -136,12 +189,18 @@
ds.arr([0.0, -1.0, 0.0], 'dimensionless'))
cam.set_width = ds.arr([8.0, 8.0, 8.0], 'code_length')
cam.resolution = (800, 800)
-
ms1 = MeshSource(ds, ('connect1', 'diffused'))
ms2 = MeshSource(ds, ('connect2', 'diffused'))
-
sc.add_source(ms1)
sc.add_source(ms2)
+ im = sc.render()
+ return compare(ds, im, "%s_composite_mesh_render" % engine)
+
+ at requires_ds(hex8)
+ at requires_module("pyembree")
+def test_composite_mesh_render_pyembree():
+ yield composite_mesh_render("embree")
- im = sc.render()
- yield compare(ds, im, "composite_mesh_render")
+ at requires_ds(hex8)
+def test_composite_mesh_render():
+ yield composite_mesh_render("yt")
diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/visualization/volume_rendering/utils.py
--- a/yt/visualization/volume_rendering/utils.py
+++ b/yt/visualization/volume_rendering/utils.py
@@ -29,7 +29,7 @@
kwargs = {'lens_type': params['lens_type']}
if engine == 'embree':
sampler = mesh_traversal.EmbreeMeshSampler(*args, **kwargs)
- elif engine == 'bvh':
+ elif engine == 'yt':
sampler = mesh_traversal.BVHMeshSampler(*args, **kwargs)
return sampler
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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-svn-spacepope.org/attachments/20160511/8ddd5f09/attachment.html>
More information about the yt-svn
mailing list