[yt-svn] commit/yt: ngoldbaum: Merged in atmyers/yt (pull request #2107)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Apr 20 11:22:30 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/248ccc5ef8b4/
Changeset: 248ccc5ef8b4
Branch: yt
User: ngoldbaum
Date: 2016-04-20 18:22:12+00:00
Summary: Merged in atmyers/yt (pull request #2107)
Allow MeshSources to optionally use the Cython ray-tracer instead of Embree.
Affected #: 14 files
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -1,6 +1,7 @@
cimport cython
import numpy as np
cimport numpy as np
+from yt.utilities.lib.element_mappings cimport ElementSampler
# ray data structure
cdef struct Ray:
@@ -11,6 +12,7 @@
np.float64_t t_near
np.float64_t t_far
np.int64_t elem_id
+ np.int64_t near_boundary
# axis-aligned bounding box
cdef struct BBox:
@@ -30,7 +32,6 @@
np.float64_t p0[3]
np.float64_t p1[3]
np.float64_t p2[3]
- np.float64_t d0, d1, d2
np.int64_t elem_id
np.float64_t centroid[3]
BBox bbox
@@ -38,6 +39,14 @@
cdef class BVH:
cdef BVHNode* root
cdef Triangle* triangles
+ 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_elem
+ cdef np.int64_t num_verts_per_elem
+ cdef np.int64_t num_field_per_elem
+ cdef ElementSampler sampler
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 intersect(self, Ray* ray) nogil
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -1,10 +1,19 @@
-cimport cython
+cimport cython
import numpy as np
cimport numpy as np
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.element_mappings cimport \
+ ElementSampler, \
+ Q1Sampler3D, \
+ P1Sampler3D, \
+ W1Sampler3D
+
+cdef ElementSampler Q1Sampler = Q1Sampler3D()
+cdef ElementSampler P1Sampler = P1Sampler3D()
+cdef ElementSampler W1Sampler = W1Sampler3D()
cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
@@ -13,7 +22,16 @@
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
@@ -60,7 +78,6 @@
if(t > DETERMINANT_EPS and t < ray.t_far):
ray.t_far = t
- ray.data_val = (1.0 - u - v)*tri.d0 + u*tri.d1 + v*tri.d2
ray.elem_id = tri.elem_id
return True
@@ -104,31 +121,62 @@
@cython.wraparound(False)
@cython.cdivision(True)
def __cinit__(self,
- np.float64_t[:, ::1] vertices,
- np.int64_t[:, ::1] indices,
- np.float64_t[:, ::1] field_data):
+ np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices,
+ np.float64_t[:, :] field_data):
- cdef np.int64_t num_elem = indices.shape[0]
- cdef np.int64_t num_tri = 12*num_elem
+ self.num_elem = indices.shape[0]
+ self.num_verts_per_elem = indices.shape[1]
+ 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.sampler = Q1Sampler
+ elif self.num_verts_per_elem == 6:
+ self.num_tri_per_elem = WEDGE_NT
+ 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.sampler = P1Sampler
+ self.num_tri = self.num_tri_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))
+
+ # create data buffers
+ cdef np.int64_t i, 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):
+ vertex_offset = i*self.num_verts_per_elem*3 + j*3
+ for k in range(3):
+ 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]
# fill our array of triangles
- cdef np.int64_t i, j, k
cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri
- self.triangles = <Triangle*> malloc(num_tri * sizeof(Triangle))
- for i in range(num_elem):
- offset = 12*i
- for j in range(12):
+ self.triangles = <Triangle*> malloc(self.num_tri * sizeof(Triangle))
+ for i in range(self.num_elem):
+ offset = self.num_tri_per_elem*i
+ for j in range(self.num_tri_per_elem):
tri_index = offset + j
tri = &(self.triangles[tri_index])
tri.elem_id = i
- v0 = indices[i][triangulate_hex[j][0]]
- v1 = indices[i][triangulate_hex[j][1]]
- v2 = indices[i][triangulate_hex[j][2]]
- tri.d0 = field_data[i][triangulate_hex[j][0]]
- tri.d1 = field_data[i][triangulate_hex[j][1]]
- tri.d2 = field_data[i][triangulate_hex[j][2]]
+ v0 = indices[i][tri_array[j][0]]
+ v1 = indices[i][tri_array[j][1]]
+ v2 = indices[i][tri_array[j][2]]
for k in range(3):
tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
@@ -137,7 +185,7 @@
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, num_tri)
+ self.root = self._recursive_build(0, self.num_tri)
cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:
@@ -148,6 +196,8 @@
def __dealloc__(self):
self._recursive_free(self.root)
free(self.triangles)
+ free(self.field_data)
+ free(self.vertices)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -189,6 +239,28 @@
@cython.cdivision(True)
cdef void intersect(self, Ray* ray) nogil:
self._recursive_intersect(ray, self.root)
+
+ if ray.elem_id < 0:
+ return
+
+ cdef np.float64_t[3] position
+ 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
+ field_ptr = self.field_data + ray.elem_id*self.num_field_per_elem
+
+ cdef np.float64_t[4] mapped_coord
+ self.sampler.map_real_to_unit(mapped_coord, vertex_ptr, position)
+ if self.num_field_per_elem == 1:
+ ray.data_val = field_ptr[0]
+ else:
+ ray.data_val = self.sampler.sample_at_unit_point(mapped_coord,
+ field_ptr)
+ ray.near_boundary = self.sampler.check_mesh_lines(mapped_coord)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -281,6 +353,7 @@
ray.t_far = INF
ray.t_near = 0.0
ray.data_val = 0
+ ray.elem_id = -1
bvh.intersect(ray)
image[i] = ray.data_val
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -35,6 +35,8 @@
double tolerance,
int direction) nogil
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil
+
cdef class P1Sampler2D(ElementSampler):
@@ -75,6 +77,8 @@
double tolerance,
int direction) nogil
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil
+
# This typedef defines a function pointer that defines the system
# of equations that will be solved by the NonlinearSolveSamplers.
@@ -170,6 +174,8 @@
double tolerance,
int direction) nogil
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil
+
cdef class W1Sampler3D(NonlinearSolveSampler3D):
@@ -190,6 +196,9 @@
double tolerance,
int direction) nogil
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil
+
+
cdef class S2Sampler3D(NonlinearSolveSampler3D):
@@ -210,6 +219,9 @@
double tolerance,
int direction) nogil
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil
+
+
cdef class NonlinearSolveSampler2D(ElementSampler):
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -95,6 +95,12 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ pass
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
cdef double sample_at_real_point(self,
double* vertices,
double* field_values,
@@ -265,6 +271,37 @@
return 1
return 0
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ cdef double u, v, w
+ cdef double thresh = 2.0e-2
+ if mapped_coord[0] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ w = mapped_coord[3]
+ elif mapped_coord[1] == 0:
+ u = mapped_coord[2]
+ v = mapped_coord[3]
+ w = mapped_coord[0]
+ elif mapped_coord[2] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[3]
+ w = mapped_coord[0]
+ else:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ w = mapped_coord[0]
+ if ((u < thresh) or
+ (v < thresh) or
+ (w < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh) or
+ (fabs(w - 1) < thresh)):
+ return 1
+ return -1
+
cdef class NonlinearSolveSampler3D(ElementSampler):
@@ -373,7 +410,6 @@
return 0
return 1
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -386,6 +422,23 @@
else:
return 0
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
+ return 1
+ elif (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
+ return 1
+ elif (fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
+ return 1
+ else:
+ return -1
+
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -528,6 +581,22 @@
else:
return 0
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
+ return 1
+ elif (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
+ return 1
+ elif (fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1 and
+ fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
+ return 1
+ else:
+ return -1
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -709,7 +778,6 @@
return 0
return 1
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -726,6 +794,32 @@
return 1
return 0
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ cdef double r, s, t
+ cdef double thresh = 5.0e-2
+ r = mapped_coord[0]
+ s = mapped_coord[1]
+ t = mapped_coord[2]
+
+ cdef int near_edge_r, near_edge_s, near_edge_t
+ near_edge_r = (r < thresh) or (fabs(r + s - 1.0) < thresh)
+ near_edge_s = (s < thresh)
+ near_edge_t = fabs(fabs(mapped_coord[2]) - 1.0) < thresh
+
+ # we use ray.instID to pass back whether the ray is near the
+ # element boundary or not (used to annotate mesh lines)
+ if (near_edge_r and near_edge_s):
+ return 1
+ elif (near_edge_r and near_edge_t):
+ return 1
+ elif (near_edge_s and near_edge_t):
+ return 1
+ else:
+ return -1
+
@cython.boundscheck(False)
@cython.wraparound(False)
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/grid_traversal.pxd
--- a/yt/utilities/lib/grid_traversal.pxd
+++ b/yt/utilities/lib/grid_traversal.pxd
@@ -25,6 +25,8 @@
np.float64_t *center
np.float64_t[:,:,:] image
np.float64_t[:,:] zbuffer
+ np.int64_t[:,:] image_used
+ np.int64_t[:,:] mesh_lines
np.float64_t pdx, pdy
np.float64_t bounds[4]
np.float64_t[:,:] camera_data # position, width, unit_vec[0,2]
@@ -59,6 +61,8 @@
cdef sampler_function *sampler
cdef public object acenter, aimage, ax_vec, ay_vec
cdef public object azbuffer
+ cdef public object aimage_used
+ cdef public object amesh_lines
cdef void *supp_data
cdef np.float64_t width[3]
cdef public object lens_type
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -383,6 +383,8 @@
*args, **kwargs):
self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
cdef np.float64_t[:,:] zbuffer
+ cdef np.int64_t[:,:] image_used
+ cdef np.int64_t[:,:] mesh_lines
cdef np.float64_t[:,:] camera_data
cdef int i
@@ -393,6 +395,10 @@
zbuffer = kwargs.pop("zbuffer", None)
if zbuffer is None:
zbuffer = np.ones((image.shape[0], image.shape[1]), "float64")
+
+ image_used = np.zeros((image.shape[0], image.shape[1]), "int64")
+ mesh_lines = np.zeros((image.shape[0], image.shape[1]), "int64")
+
self.lens_type = kwargs.pop("lens_type", None)
if self.lens_type == "plane-parallel":
self.extent_function = calculate_extent_plane_parallel
@@ -400,7 +406,7 @@
else:
if not (vp_pos.shape[0] == vp_dir.shape[0] == image.shape[0]) or \
not (vp_pos.shape[1] == vp_dir.shape[1] == image.shape[1]):
- msg = "Bad lense shape / direction for %s\n" % (self.lens_type)
+ msg = "Bad lens shape / direction for %s\n" % (self.lens_type)
msg += "Shapes: (%s - %s - %s) and (%s - %s - %s)" % (
vp_pos.shape[0], vp_dir.shape[0], image.shape[0],
vp_pos.shape[1], vp_dir.shape[1], image.shape[1])
@@ -426,7 +432,9 @@
self.image.x_vec = <np.float64_t *> x_vec.data
self.ay_vec = y_vec
self.image.y_vec = <np.float64_t *> y_vec.data
- self.image.zbuffer = zbuffer
+ self.image.zbuffer = self.azbuffer = zbuffer
+ self.image.image_used = self.aimage_used = image_used
+ self.image.mesh_lines = self.amesh_lines = mesh_lines
self.image.nv[0] = image.shape[0]
self.image.nv[1] = image.shape[1]
for i in range(4): self.image.bounds[i] = bounds[i]
@@ -501,6 +509,7 @@
self.image.vp_dir = None
self.image.zbuffer = None
self.image.camera_data = None
+ self.image.image_used = None
free(self.image)
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/mesh_construction.pxd
--- a/yt/utilities/lib/mesh_construction.pxd
+++ b/yt/utilities/lib/mesh_construction.pxd
@@ -11,6 +11,7 @@
int* element_indices # which vertices belong to which *element*
int tpe # the number of triangles per element
int vpe # the number of vertices per element
+ int fpe # the number of field values per element
ctypedef struct Patch:
float[8][3] v
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -25,7 +25,6 @@
from mesh_samplers cimport \
sample_hex, \
sample_tetra, \
- sample_element, \
sample_wedge
from pyembree.rtcore cimport \
Vertex, \
@@ -104,7 +103,7 @@
np.ndarray indices,
np.ndarray data):
- # We need now to figure out what kind of elements we've been handed.
+ # We need to figure out what kind of elements we've been handed.
if indices.shape[1] == 8:
self.vpe = HEX_NV
self.tpe = HEX_NT
@@ -186,19 +185,18 @@
datac.element_indices = self.element_indices
datac.tpe = self.tpe
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):
- if self.fpe == 1:
- self.filter_func = <rtcg.RTCFilterFunc> sample_element
- elif self.fpe == 4:
+ if self.vpe == 4:
self.filter_func = <rtcg.RTCFilterFunc> sample_tetra
- elif self.fpe == 6:
+ elif self.vpe == 6:
self.filter_func = <rtcg.RTCFilterFunc> sample_wedge
- elif self.fpe == 8:
+ elif self.vpe == 8:
self.filter_func = <rtcg.RTCFilterFunc> sample_hex
else:
raise NotImplementedError("Sampler type not implemented.")
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/mesh_samplers.pxd
--- a/yt/utilities/lib/mesh_samplers.pxd
+++ b/yt/utilities/lib/mesh_samplers.pxd
@@ -13,6 +13,3 @@
cdef void sample_hex20(void* userPtr,
rtcr.RTCRay& ray) nogil
-
-cdef void sample_element(void* userPtr,
- rtcr.RTCRay& ray) nogil
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/mesh_samplers.pyx
--- a/yt/utilities/lib/mesh_samplers.pyx
+++ b/yt/utilities/lib/mesh_samplers.pyx
@@ -97,7 +97,9 @@
for i in range(8):
element_indices[i] = data.element_indices[elem_id*8+i]
- field_data[i] = data.field_data[elem_id*8+i]
+
+ for i in range(data.fpe):
+ field_data[i] = data.field_data[elem_id*data.fpe+i]
for i in range(8):
vertices[i*3] = data.vertices[element_indices[i]].x
@@ -107,22 +109,16 @@
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
Q1Sampler.map_real_to_unit(mapped_coord, vertices, position)
- val = Q1Sampler.sample_at_unit_point(mapped_coord, field_data)
+ if data.fpe == 1:
+ val = field_data[0]
+ else:
+ val = Q1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
# we use ray.instID to pass back whether the ray is near the
# element boundary or not (used to annotate mesh lines)
- if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
- ray.instID = 1
- elif (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
- ray.instID = 1
- elif (fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
- ray.instID = 1
- else:
- ray.instID = -1
+ ray.instID = Q1Sampler.check_mesh_lines(mapped_coord)
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -152,7 +148,9 @@
for i in range(6):
element_indices[i] = data.element_indices[elem_id*6+i]
- field_data[i] = data.field_data[elem_id*6+i]
+
+ for i in range(data.fpe):
+ field_data[i] = data.field_data[elem_id*data.fpe+i]
for i in range(6):
vertices[i*3] = data.vertices[element_indices[i]].x
@@ -162,31 +160,12 @@
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
W1Sampler.map_real_to_unit(mapped_coord, vertices, position)
- val = W1Sampler.sample_at_unit_point(mapped_coord, field_data)
+ if data.fpe == 1:
+ val = field_data[0]
+ else:
+ val = W1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
-
- cdef double r, s, t
- cdef double thresh = 5.0e-2
- r = mapped_coord[0]
- s = mapped_coord[1]
- t = mapped_coord[2]
-
- cdef int near_edge_r, near_edge_s, near_edge_t
- near_edge_r = (r < thresh) or (fabs(r + s - 1.0) < thresh)
- near_edge_s = (s < thresh)
- near_edge_t = fabs(fabs(mapped_coord[2]) - 1.0) < thresh
-
- # we use ray.instID to pass back whether the ray is near the
- # element boundary or not (used to annotate mesh lines)
- if (near_edge_r and near_edge_s):
- ray.instID = 1
- elif (near_edge_r and near_edge_t):
- ray.instID = 1
- elif (near_edge_s and near_edge_t):
- ray.instID = 1
- else:
- ray.instID = -1
-
+ ray.instID = W1Sampler.check_mesh_lines(mapped_coord)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -222,21 +201,8 @@
S2Sampler.map_real_to_unit(mapped_coord, patch.vertices, position)
val = S2Sampler.sample_at_unit_point(mapped_coord, patch.field_data)
ray.time = val
-
- # we use ray.instID to pass back whether the ray is near the
- # element boundary or not (used to annotate mesh lines)
- if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
- ray.instID = 1
- elif (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
- ray.instID = 1
- elif (fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1 and
- fabs(fabs(mapped_coord[2]) - 1.0) < 1e-1):
- ray.instID = 1
- else:
- ray.instID = -1
-
+ ray.instID = S2Sampler.check_mesh_lines(mapped_coord)
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -267,54 +233,19 @@
for i in range(4):
element_indices[i] = data.element_indices[elem_id*4+i]
- field_data[i] = data.field_data[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
+ for i in range(data.fpe):
+ field_data[i] = data.field_data[elem_id*data.fpe+i]
+
# we use ray.time to pass the value of the field
cdef double mapped_coord[4]
P1Sampler.map_real_to_unit(mapped_coord, vertices, position)
- val = P1Sampler.sample_at_unit_point(mapped_coord, field_data)
+ if data.fpe == 1:
+ val = field_data[0]
+ else:
+ val = P1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
-
- cdef double u, v, w
- cdef double thresh = 2.0e-2
- u = ray.u
- v = ray.v
- w = 1.0 - u - v
- # we use ray.instID to pass back whether the ray is near the
- # element boundary or not (used to annotate mesh lines)
- if ((u < thresh) or
- (v < thresh) or
- (w < thresh) or
- (fabs(u - 1) < thresh) or
- (fabs(v - 1) < thresh) or
- (fabs(w - 1) < thresh)):
- ray.instID = 1
- else:
- ray.instID = -1
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void sample_element(void* userPtr,
- rtcr.RTCRay& ray) nogil:
- cdef int ray_id, elem_id
- cdef double val
- cdef MeshDataContainer* data
-
- data = <MeshDataContainer*> userPtr
- ray_id = ray.primID
- if ray_id == -1:
- return
-
- # 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.
- elem_id = ray_id / data.tpe
-
- val = data.field_data[elem_id]
- ray.time = val
+ ray.instID = P1Sampler.check_mesh_lines(mapped_coord)
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/mesh_traversal.pyx
--- a/yt/utilities/lib/mesh_traversal.pyx
+++ b/yt/utilities/lib/mesh_traversal.pyx
@@ -1,6 +1,6 @@
"""
-This file contains the MeshSampler class, which handles casting rays at a
-MeshSource using pyembree.
+This file contains the MeshSampler classes, which handles casting rays at a
+mesh source using either pyembree or the cython ray caster.
"""
@@ -25,6 +25,7 @@
ImageContainer
from cython.parallel import prange, parallel, threadid
from yt.visualization.image_writer import apply_colormap
+from yt.utilities.lib.bounding_volume_hierarchy cimport BVH, Ray
rtc.rtcInit(NULL)
rtc.rtcSetErrorFunction(error_printer)
@@ -42,11 +43,8 @@
def __dealloc__(self):
rtcs.rtcDeleteScene(self.scene_i)
-cdef class MeshSampler(ImageSampler):
- cdef public object image_used
- cdef public object mesh_lines
- cdef public object zbuffer
+cdef class EmbreeMeshSampler(ImageSampler):
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -70,15 +68,10 @@
cdef np.float64_t width[3]
for i in range(3):
width[i] = self.width[i]
- cdef np.ndarray[np.float64_t, ndim=1] data
cdef np.ndarray[np.int64_t, ndim=1] used
nx = im.nv[0]
ny = im.nv[1]
size = nx * ny
- used = np.empty(size, dtype="int64")
- mesh = np.empty(size, dtype="int64")
- data = np.empty(size, dtype="float64")
- zbuffer = np.empty(size, dtype="float64")
cdef rtcr.RTCRay ray
v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
@@ -99,13 +92,65 @@
ray.time = 0
ray.Ng[0] = 1e37 # we use this to track the hit distance
rtcs.rtcIntersect(scene.scene_i, ray)
- data[j] = ray.time
- used[j] = ray.primID
- mesh[j] = ray.instID
- zbuffer[j] = ray.tfar
- self.aimage = data
- self.image_used = used
- self.mesh_lines = mesh
- self.zbuffer = zbuffer
+ im.image[vi, vj, 0] = ray.time
+ im.image_used[vi, vj] = ray.primID
+ im.mesh_lines[vi, vj] = ray.instID
+ im.zbuffer[vi, vj] = ray.tfar
free(v_pos)
free(v_dir)
+
+
+cdef class BVHMeshSampler(ImageSampler):
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __call__(self,
+ BVH bvh,
+ int num_threads = 0):
+ '''
+
+ This function is supposed to cast the rays and return the
+ image.
+
+ '''
+
+ cdef int vi, vj, i, j
+ cdef ImageContainer *im = self.image
+ cdef np.float64_t *v_pos
+ cdef np.float64_t *v_dir
+ cdef np.int64_t nx, ny, size
+ cdef np.float64_t width[3]
+ for i in range(3):
+ width[i] = self.width[i]
+ cdef np.ndarray[np.float64_t, ndim=1] data
+ cdef np.ndarray[np.int64_t, ndim=1] used
+ nx = im.nv[0]
+ ny = im.nv[1]
+ size = nx * ny
+ cdef Ray* ray
+ with nogil, parallel():
+ ray = <Ray *> malloc(sizeof(Ray))
+ v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ for j in prange(size):
+ vj = j % ny
+ vi = (j - vj) / ny
+ vj = vj
+ self.vector_function(im, vi, vj, width, v_dir, v_pos)
+ for i in range(3):
+ ray.origin[i] = v_pos[i]
+ ray.direction[i] = v_dir[i]
+ ray.inv_dir[i] = 1.0 / v_dir[i]
+ ray.t_far = 1e37
+ ray.t_near = 0.0
+ ray.data_val = 0
+ ray.elem_id = -1
+ bvh.intersect(ray)
+ im.image[vi, vj, 0] = ray.data_val
+ im.image_used[vi, vj] = ray.elem_id
+ im.mesh_lines[vi, vj] = ray.near_boundary
+ im.zbuffer[vi, vj] = ray.t_far
+ free(v_pos)
+ free(v_dir)
+ free(ray)
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -20,6 +20,7 @@
from yt.utilities.lib.fp_utils cimport fmin, fmax, i64min, i64max, imin, imax, fabs
from yt.utilities.exceptions import YTPixelizeError, \
YTElementTypeNotRecognized
+from libc.stdlib cimport malloc, free
from vec3_ops cimport dot, cross, subtract
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
@@ -30,9 +31,6 @@
Q1Sampler2D, \
W1Sampler3D
-cdef extern from "stdlib.h":
- # NOTE that size_t might not be int
- void *alloca(int)
cdef extern from "pixelization_constants.h":
enum:
@@ -573,8 +571,7 @@
cdef int nvertices = conn.shape[1]
cdef int ndim = coords.shape[1]
cdef int num_field_vals = field.shape[1]
- cdef double* mapped_coord
- cdef int num_mapped_coords
+ cdef double[4] mapped_coord
cdef ElementSampler sampler
# Pick the right sampler and allocate storage for the mapped coordinate
@@ -617,10 +614,8 @@
raise RuntimeError
# allocate temporary storage
- num_mapped_coords = sampler.num_mapped_coords
- mapped_coord = <double*> alloca(sizeof(double) * num_mapped_coords)
- vertices = <np.float64_t *> alloca(ndim * sizeof(np.float64_t) * nvertices)
- field_vals = <np.float64_t *> alloca(sizeof(np.float64_t) * num_field_vals)
+ vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
+ field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
# fill the image bounds and pixel size informaton here
for i in range(ndim):
@@ -691,4 +686,7 @@
img[pi, pj, pk] = 1.0
elif sampler.check_near_edge(mapped_coord, thresh, yax):
img[pi, pj, pk] = 1.0
+
+ free(vertices)
+ free(field_vals)
return img
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -24,6 +24,7 @@
from .utils import new_volume_render_sampler, data_source_or_all, \
get_corners, new_projection_sampler, new_mesh_sampler, \
new_interpolated_projection_sampler
+from yt.utilities.lib.bounding_volume_hierarchy import BVH
from yt.visualization.image_writer import apply_colormap
from yt.data_objects.image_array import ImageArray
from .zbuffer_array import ZBuffer
@@ -353,8 +354,9 @@
self.data_source = data_source_or_all(data_source)
field = self.data_source._determine_fields(field)[0]
self.field = field
- self.mesh = None
+ self.volume = None
self.current_image = None
+ self.engine = 'embree'
# default color map
self._cmap = ytcfg.get("yt", "default_colormap")
@@ -369,8 +371,14 @@
assert(self.field is not None)
assert(self.data_source is not None)
- self.scene = mesh_traversal.YTEmbreeScene()
- self.build_mesh()
+ if self.engine == 'embree':
+ self.volume = mesh_traversal.YTEmbreeScene()
+ self.build_volume_embree()
+ elif self.engine == 'bvh':
+ self.build_volume_bvh()
+ else:
+ raise NotImplementedError("Invalid ray-tracing engine selected. "
+ "Choices are 'embree' and 'bvh'.")
def cmap():
'''
@@ -410,12 +418,60 @@
def _validate(self):
"""Make sure that all dependencies have been met"""
if self.data_source is None:
- raise RuntimeError("Data source not initialized")
+ raise RuntimeError("Data source not initialized.")
- if self.mesh is None:
- raise RuntimeError("Mesh not initialized")
+ if self.volume is None:
+ raise RuntimeError("Volume not initialized.")
- def build_mesh(self):
+ def build_volume_embree(self):
+ """
+
+ This constructs the mesh that will be ray-traced by pyembree.
+
+ """
+ ftype, fname = self.field
+ mesh_id = int(ftype[-1]) - 1
+ index = self.data_source.ds.index
+ offset = index.meshes[mesh_id]._index_offset
+ field_data = self.data_source[self.field].d # strip units
+
+ vertices = index.meshes[mesh_id].connectivity_coords
+ indices = index.meshes[mesh_id].connectivity_indices - offset
+
+ # if this is an element field, promote to 2D here
+ 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
+ # low-order geometry. Right now, high-order geometry is only
+ # implemented for 20-point hexes.
+ if indices.shape[1] == 20:
+ self.mesh = mesh_construction.QuadraticElementMesh(self.volume,
+ vertices,
+ indices,
+ field_data)
+ else:
+ # if this is another type of higher-order element, we demote
+ # to 1st order here, for now.
+ if indices.shape[1] == 27:
+ # hexahedral
+ mylog.warning("27-node hexes not yet supported, " +
+ "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,
+ indices,
+ field_data)
+
+ def build_volume_bvh(self):
"""
This constructs the mesh that will be ray-traced.
@@ -436,32 +492,27 @@
# 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.
+ # supported by Embree.
if indices.shape[1] == 20:
- self.mesh = mesh_construction.QuadraticElementMesh(self.scene,
- vertices,
- indices,
- field_data)
- else:
- # if this is another type of higher-order element, we demote
- # to 1st order here, for now.
- if indices.shape[1] == 27:
- # hexahedral
- mylog.warning("27-node hexes not yet supported, " +
- "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]
+ # 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:
+ # hexahedral
+ mylog.warning("27-node hexes not yet supported, " +
+ "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.scene,
- vertices,
- indices,
- field_data)
+ self.volume = BVH(vertices, indices, field_data)
def render(self, camera, zbuffer=None):
"""Renders an image using the provided camera
@@ -495,18 +546,17 @@
zbuffer.z.reshape(shape[:2]))
self.zbuffer = zbuffer
- self.sampler = new_mesh_sampler(camera, self)
+ self.sampler = new_mesh_sampler(camera, self, engine=self.engine)
mylog.debug("Casting rays")
- self.sampler(self.scene)
+ self.sampler(self.volume)
mylog.debug("Done casting rays")
self.finalize_image(camera)
- self.data = self.sampler.aimage
self.current_image = self.apply_colormap()
zbuffer += ZBuffer(self.current_image.astype('float64'),
- self.sampler.zbuffer)
+ self.sampler.azbuffer)
zbuffer.rgba = ImageArray(zbuffer.rgba)
self.zbuffer = zbuffer
self.current_image = self.zbuffer.rgba
@@ -523,16 +573,13 @@
# reshape data
Nx = camera.resolution[0]
Ny = camera.resolution[1]
- sam.aimage = sam.aimage.reshape(Nx, Ny)
- sam.image_used = sam.image_used.reshape(Nx, Ny)
- sam.mesh_lines = sam.mesh_lines.reshape(Nx, Ny)
- sam.zbuffer = sam.zbuffer.reshape(Nx, Ny)
+ self.data = sam.aimage[:,:,0].reshape(Nx, Ny)
# rotate
- sam.aimage = np.rot90(sam.aimage, k=2)
- sam.image_used = np.rot90(sam.image_used, k=2)
- sam.mesh_lines = np.rot90(sam.mesh_lines, k=2)
- sam.zbuffer = np.rot90(sam.zbuffer, k=2)
+ self.data = np.rot90(self.data, k=2)
+ sam.aimage_used = np.rot90(sam.aimage_used, k=2)
+ sam.amesh_lines = np.rot90(sam.amesh_lines, k=2)
+ sam.azbuffer = np.rot90(sam.azbuffer, k=2)
def annotate_mesh_lines(self, color=None, alpha=1.0):
r"""
@@ -558,7 +605,7 @@
if color is None:
color = np.array([0, 0, 0, alpha])
- locs = [self.sampler.mesh_lines == 1]
+ locs = [self.sampler.amesh_lines == 1]
self.current_image[:, :, 0][locs] = color[0]
self.current_image[:, :, 1][locs] = color[1]
@@ -592,7 +639,7 @@
color_bounds=self._color_bounds,
cmap_name=self._cmap)/255.
alpha = image[:, :, 3]
- alpha[self.sampler.image_used == -1] = 0.0
+ alpha[self.sampler.aimage_used == -1] = 0.0
image[:, :, 3] = alpha
return image
diff -r 5d118c06fa7eecc494eff5d231ebf3c47d9d4ab9 -r 248ccc5ef8b46d83682344d249c8de1c964759ab yt/visualization/volume_rendering/utils.py
--- a/yt/visualization/volume_rendering/utils.py
+++ b/yt/visualization/volume_rendering/utils.py
@@ -14,7 +14,7 @@
return data_source
-def new_mesh_sampler(camera, render_source):
+def new_mesh_sampler(camera, render_source, engine):
params = ensure_code_unit_params(camera._get_sampler_params(render_source))
args = (
np.atleast_3d(params['vp_pos']),
@@ -27,7 +27,10 @@
params['width'],
)
kwargs = {'lens_type': params['lens_type']}
- sampler = mesh_traversal.MeshSampler(*args, **kwargs)
+ if engine == 'embree':
+ sampler = mesh_traversal.EmbreeMeshSampler(*args, **kwargs)
+ elif engine == 'bvh':
+ 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/20160420/212ffecc/attachment-0001.htm>
More information about the yt-svn
mailing list