<html><body>
<p>22 new commits in yt:</p>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/bee6ed6f0887/">https://bitbucket.org/yt_analysis/yt/commits/bee6ed6f0887/</a> Changeset: bee6ed6f0887 Branch: yt User: atmyers Date: 2016-04-22 18:01:42+00:00 Summary: Fix the mesh render unit test to clear the Scene each loop iteration. Affected #: 2 files</p>
<p>diff -r 843a342ee5101c18f9acf8a490bdea39d4ee0fc2 -r bee6ed6f088779480160facd5b9a6de82b3b42c3 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -249,7 +249,7 @@</p>
<pre>position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
cdef np.float64_t* vertex_ptr</pre>
<ul><li><p>cdef np.float64_t* field_ptr</p></li></ul>
<p>+ cdef np.float64_t* field_ptr</p>
<pre> 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
</pre>
<p>diff -r 843a342ee5101c18f9acf8a490bdea39d4ee0fc2 -r bee6ed6f088779480160facd5b9a6de82b3b42c3 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 @@ -31,8 +31,8 @@</p>
<pre> images = []
ds = fake_tetrahedral_ds()</pre>
<ul><li><p>sc = Scene() for field in ds.field_list:</p></li></ul>
<p>+ sc = Scene()</p>
<pre>sc.add_source(MeshSource(ds, field))
sc.add_camera()
im = sc.render()</pre>
<p>@@ -40,6 +40,7 @@</p>
<pre>ds = fake_hexahedral_ds()
for field in ds.field_list:</pre>
<p>+ sc = Scene()</p>
<pre>sc.add_source(MeshSource(ds, field))
sc.add_camera()
im = sc.render()</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/68e766e928a3/">https://bitbucket.org/yt_analysis/yt/commits/68e766e928a3/</a> Changeset: 68e766e928a3 Branch: yt User: atmyers Date: 2016-04-22 21:07:21+00:00 Summary: missing a dep in setup.py Affected #: 1 file</p>
<p>diff -r bee6ed6f088779480160facd5b9a6de82b3b42c3 -r 68e766e928a3328b90c8512156b67e438d5030ad setup.py --- a/setup.py +++ b/setup.py @@ -279,7 +279,8 @@</p>
<pre> Extension("yt.utilities.lib.mesh_traversal",
["yt/utilities/lib/mesh_traversal.pyx"],
depends=["yt/utilities/lib/mesh_traversal.pxd",</pre>
<ul><li><p>"yt/utilities/lib/grid_traversal.pxd"]),</p></li></ul>
<p>+ “yt/utilities/lib/grid_traversal.pxd”, + "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),</p>
<pre> Extension("yt.utilities.lib.mesh_samplers",
["yt/utilities/lib/mesh_samplers.pyx"],
depends=["yt/utilities/lib/mesh_samplers.pxd",</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/04c7d22d4aeb/">https://bitbucket.org/yt_analysis/yt/commits/04c7d22d4aeb/</a> Changeset: 04c7d22d4aeb Branch: yt User: atmyers Date: 2016-04-22 21:49:22+00:00 Summary: seperate out triangle centroid and triangle bounding box operations into separate functions Affected #: 1 file</p>
<p>diff -r 68e766e928a3328b90c8512156b67e438d5030ad -r 04c7d22d4aeb14cd7df92ae3f07d5c5acb785be7 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -42,9 +42,13 @@</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef np.int64_t ray_triangle_intersect(Ray* ray, const Triangle* tri) nogil: +cdef np.int64_t ray_triangle_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil:</p>
<pre># https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
</pre>
<p>+ cdef Triangle tri = (<Triangle*> primitives)[item] +</p>
<pre># edge vectors
cdef np.float64_t e1[3]
cdef np.float64_t e2[3]</pre>
<p>@@ -87,6 +91,33 @@</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>+cdef void triangle_centroid(const void <strong>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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void triangle_bbox(const void *primitives, + const np.int64_t item, + BBox</strong> 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]) + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True)</p>
<pre>cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil:
# https://tavianator.com/fast-branchless-raybounding-box-intersections/
</pre>
<p>@@ -181,9 +212,8 @@</p>
<pre>tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>tri.centroid[k] = (tri.p0[k] + tri.p1[k] + tri.p2[k]) / 3.0</p></li>
<li><p>tri.bbox.left_edge[k] = fmin(fmin(tri.p0[k], tri.p1[k]), tri.p2[k])</p></li>
<li><p>tri.bbox.right_edge[k] = fmax(fmax(tri.p0[k], tri.p1[k]), tri.p2[k])</p></li></ul>
<p>+ triangle_centroid(self.triangles, tri_index, tri.centroid) + triangle_bbox(self.triangles, tri_index, &(tri.bbox))</p>
<pre> self.root = self._recursive_build(0, self.num_tri)
</pre>
<p>@@ -276,8 +306,7 @@</p>
<pre> cdef Triangle* tri
if (node.end - node.begin) <= LEAF_SIZE:
for i in range(node.begin, node.end):</pre>
<ul><li><p>tri = &(self.triangles[i])</p></li>
<li><p>hit = ray_triangle_intersect(ray, tri)</p></li></ul>
<p>+ hit = ray_triangle_intersect(self.triangles, i, ray)</p>
<pre> return
# if not leaf, intersect with left and right children</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/4b589cc79bf7/">https://bitbucket.org/yt_analysis/yt/commits/4b589cc79bf7/</a> Changeset: 4b589cc79bf7 Branch: yt User: atmyers Date: 2016-04-22 22:08:10+00:00 Summary: seperate bounding boxes and centroids from the geometry data structure Affected #: 2 files</p>
<p>diff -r 04c7d22d4aeb14cd7df92ae3f07d5c5acb785be7 -r 4b589cc79bf7537245847a7567d105f90a7b2664 yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -33,12 +33,12 @@</p>
<pre>np.float64_t p1[3]
np.float64_t p2[3]
np.int64_t elem_id</pre>
<ul><li><p>np.float64_t centroid[3]</p></li>
<li><p>BBox bbox</p></li></ul>
<pre>cdef class BVH:
cdef BVHNode* root
cdef Triangle* triangles</pre>
<p>+ cdef np.float64_t** centroids + cdef BBox* bboxes</p>
<pre>cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef np.int64_t num_tri_per_elem</pre>
<p>diff -r 04c7d22d4aeb14cd7df92ae3f07d5c5acb785be7 -r 4b589cc79bf7537245847a7567d105f90a7b2664 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -199,6 +199,10 @@</p>
<pre>cdef np.int64_t v0, v1, v2
cdef Triangle* tri
self.triangles = <Triangle*> malloc(self.num_tri * sizeof(Triangle))</pre>
<p>+ self.centroids = <np.float64_t**> malloc(self.num_tri * sizeof(np.float64_t*)) + for i in range(self.num_tri): + self.centroids[i] = <np.float64_t*> malloc(3*sizeof(np.float64_t)) + self.bboxes = <BBox*> malloc(self.num_tri * sizeof(BBox))</p>
<pre> for i in range(self.num_elem):
offset = self.num_tri_per_elem*i
for j in range(self.num_tri_per_elem):</pre>
<p>@@ -212,8 +216,8 @@</p>
<pre>tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>triangle_centroid(self.triangles, tri_index, tri.centroid)</p></li>
<li><p>triangle_bbox(self.triangles, tri_index, &(tri.bbox))</p></li></ul>
<p>+ triangle_centroid(self.triangles, tri_index, self.centroids[tri_index]) + triangle_bbox(self.triangles, tri_index, &(self.bboxes[tri_index]))</p>
<pre> self.root = self._recursive_build(0, self.num_tri)
</pre>
<p>@@ -226,6 +230,11 @@</p>
<pre>def __dealloc__(self):
self._recursive_free(self.root)
free(self.triangles)</pre>
<p>+ cdef np.int64_t i + for i in range(self.num_tri): + free(self.centroids[i]) + free(self.centroids) + free(self.bboxes)</p>
<pre> free(self.field_data)
free(self.vertices)
</pre>
<p>@@ -240,11 +249,15 @@</p>
<pre># will have centroids *greater* than "split" along "ax".
cdef np.int64_t mid = begin
while (begin != end):</pre>
<ul><li><p>if self.triangles[mid].centroid[ax] > split:</p></li></ul>
<p>+ if self.centroids[mid][ax] > split:</p>
<pre>mid += 1</pre>
<ul><li><p>elif self.triangles[begin].centroid[ax] > split:</p></li></ul>
<p>+ elif self.centroids[begin][ax] > split:</p>
<pre>self.triangles[mid], self.triangles[begin] = \
self.triangles[begin], self.triangles[mid]</pre>
<p>+ self.centroids[mid], self.centroids[begin] = \ + self.centroids[begin], self.centroids[mid] + self.bboxes[mid], self.bboxes[begin] = \ + self.bboxes[begin], self.bboxes[mid]</p>
<pre> mid += 1
begin += 1
return mid</pre>
<p>@@ -255,13 +268,13 @@</p>
<pre> cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil:
cdef np.int64_t i, j</pre>
<ul><li><p>cdef BBox box = self.triangles[begin].bbox</p></li></ul>
<p>+ cdef BBox box = self.bboxes[begin]</p>
<pre> for i in range(begin+1, end):
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],</pre>
<ul><li><p>self.triangles[i].bbox.left_edge[j])</p></li></ul>
<p>+ self.bboxes[i].left_edge[j])</p>
<pre>box.right_edge[j] = fmax(box.right_edge[j],</pre>
<ul><li><p>self.triangles[i].bbox.right_edge[j])</p></li></ul>
<p>+ self.bboxes[i].right_edge[j])</p>
<pre> node.bbox = box
@cython.boundscheck(False)</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/5466aa0c7789/">https://bitbucket.org/yt_analysis/yt/commits/5466aa0c7789/</a> Changeset: 5466aa0c7789 Branch: yt User: atmyers Date: 2016-04-22 22:35:00+00:00 Summary: putting geometric primitive functions and data structures into a separate file. Affected #: 5 files</p>
<p>diff -r 4b589cc79bf7537245847a7567d105f90a7b2664 -r 5466aa0c7789146434efca385d7395c49d2d013a setup.py --- a/setup.py +++ b/setup.py @@ -157,6 +157,12 @@</p>
<pre> Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),</pre>
<p>+ 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"]),</p>
<pre> Extension("yt.utilities.lib.bounding_volume_hierarchy",
["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
include_dirs=["yt/utilities/lib/"],</pre>
<p>@@ -164,7 +170,7 @@</p>
<pre>extra_link_args=omp_args,
libraries=std_libs,
depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",</pre>
<ul><li><p>"yt/utilities/lib/vec3_ops.pxd"]),</p></li></ul>
<p>+ "yt/utilities/lib/primitives.pxd"]),</p>
<pre> Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],
include_dirs=["yt/utilities/lib/",</pre>
<p>diff -r 4b589cc79bf7537245847a7567d105f90a7b2664 -r 5466aa0c7789146434efca385d7395c49d2d013a yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -2,6 +2,7 @@</p>
<pre>import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler</pre>
<p>+from yt.utilities.lib.primitives cimport Triangle</p>
<pre># ray data structure
cdef struct Ray:</pre>
<p>@@ -27,13 +28,6 @@</p>
<pre>BVHNode* right
BBox bbox
</pre>
<p>-# triangle data structure -cdef struct Triangle:</p>
<ul><li><p>np.float64_t p0[3]</p></li>
<li><p>np.float64_t p1[3]</p></li>
<li><p>np.float64_t p2[3]</p></li>
<li><p>np.int64_t elem_id</p></li></ul>
<p>–</p>
<pre>cdef class BVH:
cdef BVHNode* root
cdef Triangle* triangles</pre>
<p>diff -r 4b589cc79bf7537245847a7567d105f90a7b2664 -r 5466aa0c7789146434efca385d7395c49d2d013a yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -4,7 +4,11 @@</p>
<pre>from libc.math cimport fabs
from libc.stdlib cimport malloc, free
from cython.parallel import parallel, prange</pre>
<p>-from vec3_ops cimport dot, subtract, cross +from yt.utilities.lib.primitives cimport \ + Triangle, \ + ray_triangle_intersect, \ + triangle_centroid, \ + triangle_bbox</p>
<pre>from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \</pre>
<p>@@ -34,7 +38,6 @@</p>
<pre> int triangulate_wedge[MAX_NUM_TRI][3]
# define some constants</pre>
<p>-cdef np.float64_t DETERMINANT_EPS = 1.0e-10</p>
<pre>cdef np.float64_t INF = np.inf
cdef np.int64_t LEAF_SIZE = 16
</pre>
<p>@@ -42,82 +45,6 @@</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef np.int64_t ray_triangle_intersect(const void* primitives,</p>
<ul><li><p>const np.int64_t item,</p></li>
<li><p>Ray* ray) nogil:</p></li></ul>
<p>-# <a href="https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm">https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm</a> –</p>
<ul><li><p>cdef Triangle tri = (<Triangle*> primitives)[item]</p></li></ul>
<p>–</p>
<ul><li><p># edge vectors</p></li>
<li><p>cdef np.float64_t e1[3]</p></li>
<li><p>cdef np.float64_t e2[3]</p></li>
<li><p>subtract(tri.p1, tri.p0, e1)</p></li>
<li><p>subtract(tri.p2, tri.p0, e2)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t P[3]</p></li>
<li><p>cross(ray.direction, e2, P)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t det, inv_det</p></li>
<li><p>det = dot(e1, P)</p></li>
<li><p>if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):</p></li>
<li><p>return False</p></li>
<li><p>inv_det = 1.0 / det</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t T[3]</p></li>
<li><p>subtract(ray.origin, tri.p0, T)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t u = dot(T, P) * inv_det</p></li>
<li><p>if(u < 0.0 or u > 1.0):</p></li>
<li><p>return False</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t Q[3]</p></li>
<li><p>cross(T, e1, Q)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t v = dot(ray.direction, Q) * inv_det</p></li>
<li><p>if(v < 0.0 or u + v > 1.0):</p></li>
<li><p>return False</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t t = dot(e2, Q) * inv_det</p></li></ul>
<p>–</p>
<ul><li><p>if(t > DETERMINANT_EPS and t < ray.t_far):</p></li>
<li><p>ray.t_far = t</p></li>
<li><p>ray.elem_id = tri.elem_id</p></li>
<li><p>return True</p></li></ul>
<p>–</p>
<ul><li><p>return False</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void triangle_centroid(const void *primitives,</p>
<ul><li><p>const np.int64_t item,</p></li>
<li><p>np.float64_t[3] centroid) nogil:</p></li></ul>
<p>–</p>
<ul><li><p>cdef Triangle tri = (<Triangle*> primitives)[item]</p></li>
<li><p>cdef np.int64_t i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>centroid[i] = (tri.p0[i] + tri.p1[i] + tri.p2[i]) / 3.0</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void triangle_bbox(const void *primitives,</p>
<ul><li><p>const np.int64_t item,</p></li>
<li><p>BBox* bbox) nogil:</p></li></ul>
<p>–</p>
<ul><li><p>cdef Triangle tri = (<Triangle*> primitives)[item]</p></li>
<li><p>cdef np.int64_t i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>bbox.left_edge[i] = fmin(fmin(tri.p0[i], tri.p1[i]), tri.p2[i])</p></li>
<li><p>bbox.right_edge[i] = fmax(fmax(tri.p0[i], tri.p1[i]), tri.p2[i])</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True)</p>
<pre>cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil:
# https://tavianator.com/fast-branchless-raybounding-box-intersections/
</pre>
<p>diff -r 4b589cc79bf7537245847a7567d105f90a7b2664 -r 5466aa0c7789146434efca385d7395c49d2d013a yt/utilities/lib/primitives.pxd --- /dev/null +++ b/yt/utilities/lib/primitives.pxd @@ -0,0 +1,25 @@ +cimport cython +import numpy as np +cimport numpy as np +from vec3_ops cimport dot, subtract, cross +from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, 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 + +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 <strong>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</strong> bbox) nogil +</p>
<p>diff -r 4b589cc79bf7537245847a7567d105f90a7b2664 -r 5466aa0c7789146434efca385d7395c49d2d013a yt/utilities/lib/primitives.pyx --- /dev/null +++ b/yt/utilities/lib/primitives.pyx @@ -0,0 +1,87 @@ +cimport cython +import numpy as np +cimport numpy as np +from vec3_ops cimport dot, subtract, cross +from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox + +cdef np.float64_t DETERMINANT_EPS = 1.0e-10 + +cdef extern from “platform_dep.h” nogil: + double fmax(double x, double y) + double fmin(double x, double y) + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline np.int64_t ray_triangle_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil: +# <a href="https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm">https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm</a> + + 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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void triangle_centroid(const void <strong>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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void triangle_bbox(const void *primitives, + const np.int64_t item, + BBox</strong> 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])</p>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/b8f6257d4bff/">https://bitbucket.org/yt_analysis/yt/commits/b8f6257d4bff/</a> Changeset: b8f6257d4bff Branch: yt User: atmyers Date: 2016-04-23 21:48:03+00:00 Summary: implementing patch geometry functions for BVH. Affected #: 2 files</p>
<p>diff -r 5466aa0c7789146434efca385d7395c49d2d013a -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 yt/utilities/lib/primitives.pxd --- a/yt/utilities/lib/primitives.pxd +++ b/yt/utilities/lib/primitives.pxd @@ -4,7 +4,6 @@</p>
<pre>from vec3_ops cimport dot, subtract, cross
from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox
</pre>
<p>-# triangle data structure</p>
<pre>cdef struct Triangle:
np.float64_t p0[3]
np.float64_t p1[3]</pre>
<p>@@ -23,3 +22,18 @@</p>
<pre> const np.int64_t item,
BBox* bbox) nogil
</pre>
<p>+cdef struct Patch: + np.float64_t[8][3] v # 8 vertices per patch + np.int64_t elem_id + +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 <strong>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</strong> bbox) nogil</p>
<p>diff -r 5466aa0c7789146434efca385d7395c49d2d013a -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -1,7 +1,9 @@</p>
<pre>cimport cython
import numpy as np
cimport numpy as np</pre>
<p>-from vec3_ops cimport dot, subtract, cross +from libc.math cimport fabs + +from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance</p>
<pre>from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox
cdef np.float64_t DETERMINANT_EPS = 1.0e-10</pre>
<p>@@ -10,7 +12,6 @@</p>
<pre> double fmax(double x, double y)
double fmin(double x, double y)
</pre>
<p>–</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>@@ -66,7 +67,7 @@</p>
<pre>cdef inline void triangle_centroid(const void *primitives,
const np.int64_t item,
np.float64_t[3] centroid) nogil:</pre>
<p>– +</p>
<pre>cdef Triangle tri = (<Triangle*> primitives)[item]
cdef np.int64_t i
for i in range(3):</pre>
<p>@@ -85,3 +86,180 @@</p>
<pre>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])</pre>
<p>+ + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void patchSurfaceFunc(const Patch* patch, + const np.float64_t u, + const np.float64_t v, + np.float64_t[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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void patchSurfaceDerivU(const Patch* patch, + const np.float64_t u, + const np.float64_t v, + np.float64_t[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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void patchSurfaceDerivV(const Patch* patch, + const np.float64_t u, + const np.float64_t v, + np.float64_t[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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline np.int64_t ray_patch_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil: + + cdef Patch patch = (<Patch*> primitives)[item] + + # first we compute the two planes that define the ray. + cdef np.float64_t[3] n, N1, N2 + cdef np.float64_t 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 np.float64_t d1 = -dot(N1, ray.origin) + cdef np.float64_t d2 = -dot(N2, ray.origin) + + # the initial guess is set to zero + cdef np.float64_t u = 0.0 + cdef np.float64_t v = 0.0 + cdef np.float64_t[3] S + patchSurfaceFunc(&patch, u, v, S) + cdef np.float64_t fu = dot(N1, S) + d1 + cdef np.float64_t fv = dot(N2, S) + d2 + cdef np.float64_t err = fmax(fabs(fu), fabs(fv)) + + # begin Newton interation + cdef np.float64_t tol = 1.0e-5 + cdef int iterations = 0 + cdef int max_iter = 10 + cdef np.float64_t[3] Su + cdef np.float64_t[3] Sv + cdef np.float64_t 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 np.float64_t t = distance(S, ray.origin) + + # only count this is it's the closest hit + if (t < ray.t_near or t > ray.t_far): + return False + + if (fabs(u) <= 1.0 and fabs(v) <= 1.0 and iterations < max_iter): + + # we have a hit, so update ray information + ray.t_far = t + ray.elem_id = patch.elem_id + return True + + return False + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void patch_centroid(const void <strong>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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void patch_bbox(const void *primitives, + const np.int64_t item, + BBox</strong> 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])</p>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/152e8577389d/">https://bitbucket.org/yt_analysis/yt/commits/152e8577389d/</a> Changeset: 152e8577389d Branch: yt User: atmyers Date: 2016-04-23 23:09:12+00:00 Summary: refactoring patch intersection code to reduce duplication. Affected #: 5 files</p>
<p>diff -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b yt/utilities/lib/mesh_intersection.pxd --- a/yt/utilities/lib/mesh_intersection.pxd +++ b/yt/utilities/lib/mesh_intersection.pxd @@ -4,15 +4,10 @@</p>
<pre>from yt.utilities.lib.mesh_construction cimport Patch
cimport cython
</pre>
<p>-cdef void patchIntersectFunc(Patch* patches,</p>
<ul><li><p>rtcr.RTCRay& ray,</p></li></ul>
<p>+cdef void patchIntersectFunc(Patch* patches, + rtcr.RTCRay& ray,</p>
<pre> size_t item) nogil
</pre>
<p>-cdef void patchBoundsFunc(Patch* patches,</p>
<ul><li><p>size_t item,</p></li></ul>
<p>+cdef void patchBoundsFunc(Patch* patches, + size_t item,</p>
<pre>rtcg.RTCBounds* bounds_o) nogil</pre>
<p>– -cdef void patchSurfaceFunc(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] S) nogil</p></li></ul>
<p>diff -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b yt/utilities/lib/mesh_intersection.pyx --- a/yt/utilities/lib/mesh_intersection.pyx +++ b/yt/utilities/lib/mesh_intersection.pyx @@ -21,66 +21,16 @@</p>
<pre>cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
from yt.utilities.lib.mesh_samplers cimport sample_hex20</pre>
<p>+from yt.utilities.lib.primitives cimport \ + patchSurfaceFunc, \ + patchSurfaceDerivU, \ + patchSurfaceDerivV</p>
<pre>from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef void patchSurfaceFunc(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] S) nogil:</p></li></ul>
<p>–</p>
<ul><li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>S[i] = 0.25*(1.0 – u)*(1.0 – v)*(-u – v – 1)*patch.v[0][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 – v)*( u – v – 1)*patch.v[1][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 + v)*( u + v – 1)*patch.v[2][i] + \</p></li>
<li><p>0.25*(1.0 – u)*(1.0 + v)*(-u + v – 1)*patch.v[3][i] + \</p></li>
<li><p>0.5*(1 – u)*(1 – v*v)*patch.v[4][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 – v)*patch.v[5][i] + \</p></li>
<li><p>0.5*(1 + u)*(1 – v*v)*patch.v[6][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 + v)*patch.v[7][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void patchSurfaceDerivU(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] Su) nogil:</p></li>
<li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>Su[i] = (-0.25*(v – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(v – 1.0)*(u – v – 1) – 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u – v + 1) + 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[4][i] + u*(v – 1.0)*patch.v[5][i] – \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[6][i] – u*(v + 1.0)*patch.v[7][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void patchSurfaceDerivV(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] Sv) nogil:</p></li>
<li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>Sv[i] = (-0.25*(u – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(u + 1.0)*(u – v – 1) + 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(u + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(u – 1.0)*(u – v + 1) – 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[5][i] + v*(u – 1.0)*patch.v[4][i] – \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[7][i] – v*(u + 1.0)*patch.v[6][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True)</p>
<pre>cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:</pre>
<p>@@ -144,7 +94,7 @@</p>
<pre>cdef float u = 0.0
cdef float v = 0.0
cdef float[3] S</pre>
<ul><li><p>patchSurfaceFunc(patch, u, v, S)</p></li></ul>
<p>+ patchSurfaceFunc(patch.v, u, v, S)</p>
<pre>cdef float fu = dot(N1, S) + d1
cdef float fv = dot(N2, S) + d2
cdef float err = fmax(fabs(fu), fabs(fv))</pre>
<p>@@ -158,8 +108,8 @@</p>
<pre>cdef float J11, J12, J21, J22, det
while ((err > tol) and (iterations < max_iter)):
# compute the Jacobian</pre>
<ul><li><p>patchSurfaceDerivU(patch, u, v, Su)</p></li>
<li><p>patchSurfaceDerivV(patch, u, v, Sv)</p></li></ul>
<p>+ patchSurfaceDerivU(patch.v, u, v, Su) + patchSurfaceDerivV(patch.v, u, v, Sv)</p>
<pre>J11 = dot(N1, Su)
J12 = dot(N1, Sv)
J21 = dot(N2, Su)</pre>
<p>@@ -170,7 +120,7 @@</p>
<pre>u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
</pre>
<ul><li><p>patchSurfaceFunc(patch, u, v, S)</p></li></ul>
<p>+ patchSurfaceFunc(patch.v, u, v, S)</p>
<pre> fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
</pre>
<p>diff -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b yt/utilities/lib/mesh_samplers.pyx --- a/yt/utilities/lib/mesh_samplers.pyx +++ b/yt/utilities/lib/mesh_samplers.pyx @@ -19,7 +19,7 @@</p>
<pre>from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
Patch</pre>
<p>-from yt.utilities.lib.mesh_intersection cimport patchSurfaceFunc +from yt.utilities.lib.primitives cimport patchSurfaceFunc</p>
<pre>from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \</pre>
<p>@@ -192,7 +192,7 @@</p>
<pre> elem_id = ray_id / 6
# fills "position" with the physical position of the hit</pre>
<ul><li><p>patchSurfaceFunc(data[ray_id], ray.u, ray.v, pos)</p></li></ul>
<p>+ patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)</p>
<pre> for i in range(3):
position[i] = <double> pos[i]
</pre>
<p>diff -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b yt/utilities/lib/primitives.pxd --- a/yt/utilities/lib/primitives.pxd +++ b/yt/utilities/lib/primitives.pxd @@ -25,6 +25,21 @@</p>
<pre>cdef struct Patch:
np.float64_t[8][3] v # 8 vertices per patch
np.int64_t elem_id</pre>
<p>+ +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</p>
<pre>cdef inline np.int64_t ray_patch_intersect(const void* primitives,
const np.int64_t item,</pre>
<p>diff -r b8f6257d4bff425d0e0a24b1e6b4599a04b55212 -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -1,6 +1,7 @@</p>
<pre>cimport cython
import numpy as np
cimport numpy as np</pre>
<p>+cimport cython.floating</p>
<pre>from libc.math cimport fabs
from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance</pre>
<p>@@ -91,55 +92,55 @@</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef void patchSurfaceFunc(const Patch* patch,</p>
<ul><li><p>const np.float64_t u,</p></li>
<li><p>const np.float64_t v,</p></li>
<li><p>np.float64_t[3] S) nogil:</p></li></ul>
<p>+cdef void patchSurfaceFunc(const cython.floating[8][3] verts, + const cython.floating u, + const cython.floating v, + cython.floating[3] S) nogil:</p>
<pre>cdef int i
for i in range(3):</pre>
<ul><li><p>S[i] = 0.25*(1.0 – u)*(1.0 – v)*(-u – v – 1)*patch.v[0][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 – v)*( u – v – 1)*patch.v[1][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 + v)*( u + v – 1)*patch.v[2][i] + \</p></li>
<li><p>0.25*(1.0 – u)*(1.0 + v)*(-u + v – 1)*patch.v[3][i] + \</p></li>
<li><p>0.5*(1 – u)*(1 – v*v)*patch.v[4][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 – v)*patch.v[5][i] + \</p></li>
<li><p>0.5*(1 + u)*(1 – v*v)*patch.v[6][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 + v)*patch.v[7][i]</p></li></ul>
<p>+ 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]</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef void patchSurfaceDerivU(const Patch* patch,</p>
<ul><li><p>const np.float64_t u,</p></li>
<li><p>const np.float64_t v,</p></li>
<li><p>np.float64_t[3] Su) nogil:</p></li></ul>
<p>+cdef void patchSurfaceDerivU(const cython.floating[8][3] verts, + const cython.floating u, + const cython.floating v, + cython.floating[3] Su) nogil:</p>
<pre>cdef int i
for i in range(3):</pre>
<ul><li><p>Su[i] = (-0.25*(v – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(v – 1.0)*(u – v – 1) – 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u – v + 1) + 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[4][i] + u*(v – 1.0)*patch.v[5][i] – \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[6][i] – u*(v + 1.0)*patch.v[7][i]</p></li></ul>
<p>+ 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]</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef void patchSurfaceDerivV(const Patch* patch,</p>
<ul><li><p>const np.float64_t u,</p></li>
<li><p>const np.float64_t v,</p></li>
<li><p>np.float64_t[3] Sv) nogil:</p></li></ul>
<p>+cdef void patchSurfaceDerivV(const cython.floating[8][3] verts, + const cython.floating u, + const cython.floating v, + cython.floating[3] Sv) nogil:</p>
<pre>cdef int i
for i in range(3):</pre>
<ul><li><p>Sv[i] = (-0.25*(u – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(u + 1.0)*(u – v – 1) + 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(u + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(u – 1.0)*(u – v + 1) – 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[5][i] + v*(u – 1.0)*patch.v[4][i] – \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[7][i] – v*(u + 1.0)*patch.v[6][i]</p></li></ul>
<p>+ 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]</p>
<pre>@cython.boundscheck(False)</pre>
<p>@@ -174,7 +175,7 @@</p>
<pre>cdef np.float64_t u = 0.0
cdef np.float64_t v = 0.0
cdef np.float64_t[3] S</pre>
<ul><li><p>patchSurfaceFunc(&patch, u, v, S)</p></li></ul>
<p>+ patchSurfaceFunc(patch.v, u, v, S)</p>
<pre>cdef np.float64_t fu = dot(N1, S) + d1
cdef np.float64_t fv = dot(N2, S) + d2
cdef np.float64_t err = fmax(fabs(fu), fabs(fv))</pre>
<p>@@ -188,8 +189,8 @@</p>
<pre>cdef np.float64_t J11, J12, J21, J22, det
while ((err > tol) and (iterations < max_iter)):
# compute the Jacobian</pre>
<ul><li><p>patchSurfaceDerivU(&patch, u, v, Su)</p></li>
<li><p>patchSurfaceDerivV(&patch, u, v, Sv)</p></li></ul>
<p>+ patchSurfaceDerivU(patch.v, u, v, Su) + patchSurfaceDerivV(patch.v, u, v, Sv)</p>
<pre>J11 = dot(N1, Su)
J12 = dot(N1, Sv)
J21 = dot(N2, Su)</pre>
<p>@@ -200,7 +201,7 @@</p>
<pre>u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
</pre>
<ul><li><p>patchSurfaceFunc(&patch, u, v, S)</p></li></ul>
<p>+ patchSurfaceFunc(patch.v, u, v, S)</p>
<pre> fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/88ad574799eb/">https://bitbucket.org/yt_analysis/yt/commits/88ad574799eb/</a> Changeset: 88ad574799eb Branch: yt User: atmyers Date: 2016-04-24 03:13:52+00:00 Summary: adding a bunch of missing dependencies to the cython extensions in setup.py. Affected #: 4 files</p>
<p>diff -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b -r 88ad574799ebf38c392490f4857d07763e9a7006 setup.py --- a/setup.py +++ b/setup.py @@ -170,6 +170,7 @@</p>
<pre>extra_link_args=omp_args,
libraries=std_libs,
depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",</pre>
<p>+ “yt/utilities/lib/element_mappings.pxd”,</p>
<pre> "yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],</pre>
<p>@@ -281,7 +282,10 @@</p>
<pre> embree_extensions = [
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],</pre>
<ul><li><p>depends=["yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ 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"]),</p>
<pre> Extension("yt.utilities.lib.mesh_traversal",
["yt/utilities/lib/mesh_traversal.pyx"],
depends=["yt/utilities/lib/mesh_traversal.pxd",</pre>
<p>@@ -291,11 +295,17 @@</p>
<pre>["yt/utilities/lib/mesh_samplers.pyx"],
depends=["yt/utilities/lib/mesh_samplers.pxd",
"yt/utilities/lib/element_mappings.pxd",</pre>
<ul><li><p>"yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ “yt/utilities/lib/mesh_construction.pxd”, + “yt/utilities/lib/bounding_volume_hierarchy.pxd”, + "yt/utilities/lib/primitives.pxd"]),</p>
<pre> Extension("yt.utilities.lib.mesh_intersection",
["yt/utilities/lib/mesh_intersection.pyx"],
depends=["yt/utilities/lib/mesh_intersection.pxd",</pre>
<ul><li><p>"yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ “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"]),</p>
<pre> ]
embree_prefix = os.path.abspath(read_embree_location())</pre>
<p>diff -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b -r 88ad574799ebf38c392490f4857d07763e9a7006 yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -27,7 +27,22 @@</p>
<pre>BVHNode* left
BVHNode* right
BBox bbox</pre>
<p>– + +# 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 <strong>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</strong> bbox) nogil +</p>
<pre>cdef class BVH:
cdef BVHNode* root
cdef Triangle* triangles</pre>
<p>@@ -41,6 +56,9 @@</p>
<pre>cdef np.int64_t num_verts_per_elem
cdef np.int64_t num_field_per_elem
cdef ElementSampler sampler</pre>
<p>+ cdef centroid_func_type get_centroid + cdef bbox_func_type get_bbox + cdef intersect_func_type get_intersect</p>
<pre> 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</pre>
<p>diff -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b -r 88ad574799ebf38c392490f4857d07763e9a7006 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -4,6 +4,7 @@</p>
<pre>from libc.math cimport fabs
from libc.stdlib cimport malloc, free
from cython.parallel import parallel, prange</pre>
<p>+</p>
<pre>from yt.utilities.lib.primitives cimport \
Triangle, \
ray_triangle_intersect, \</pre>
<p>@@ -87,6 +88,10 @@</p>
<pre> self.num_verts_per_elem = indices.shape[1]
self.num_field_per_elem = field_data.shape[1]
</pre>
<p>+ self.get_centroid = triangle_centroid + self.get_bbox = triangle_bbox + self.get_intersect = ray_triangle_intersect +</p>
<pre># 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:</pre>
<p>@@ -121,7 +126,7 @@</p>
<pre> for j in range(self.num_field_per_elem):
self.field_data[field_offset + j] = field_data[i][j]
</pre>
<ul><li><p># fill our array of triangles</p></li></ul>
<p>+ # fill our array of primitives</p>
<pre>cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri</pre>
<p>@@ -143,8 +148,12 @@</p>
<pre>tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>triangle_centroid(self.triangles, tri_index, self.centroids[tri_index])</p></li>
<li><p>triangle_bbox(self.triangles, tri_index, &(self.bboxes[tri_index]))</p></li></ul>
<p>+ self.get_centroid(self.triangles, + tri_index, + self.centroids[tri_index]) + self.get_bbox(self.triangles, + tri_index, + &(self.bboxes[tri_index]))</p>
<pre> self.root = self._recursive_build(0, self.num_tri)
</pre>
<p>@@ -171,8 +180,8 @@</p>
<pre> 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</pre>
<ul><li><p># to the left of mid have centroids less than or equal to “split”</p></li>
<li><p># along the direction “ax”. All the triangles to the right of mid</p></li></ul>
<p>+ # 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</p>
<pre># will have centroids *greater* than "split" along "ax".
cdef np.int64_t mid = begin
while (begin != end):</pre>
<p>@@ -246,7 +255,7 @@</p>
<pre> cdef Triangle* tri
if (node.end - node.begin) <= LEAF_SIZE:
for i in range(node.begin, node.end):</pre>
<ul><li><p>hit = ray_triangle_intersect(self.triangles, i, ray)</p></li></ul>
<p>+ hit = self.get_intersect(self.triangles, i, ray)</p>
<pre> return
# if not leaf, intersect with left and right children</pre>
<p>diff -r 152e8577389d8a3c4970d7b2f5f3a3970a177b9b -r 88ad574799ebf38c392490f4857d07763e9a7006 yt/utilities/lib/mesh_construction.pyx --- a/yt/utilities/lib/mesh_construction.pyx +++ b/yt/utilities/lib/mesh_construction.pyx @@ -17,25 +17,27 @@</p>
<pre>import numpy as np
cimport cython</pre>
<p>+from libc.stdlib cimport malloc, free +from libc.math cimport fmax, sqrt +cimport numpy as np +</p>
<pre>cimport pyembree.rtcore as rtc</pre>
<p>-from mesh_traversal cimport YTEmbreeScene</p>
<pre>cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu</pre>
<p>+from pyembree.rtcore cimport \ + Vertex, \ + Triangle, \ + Vec3f + +from mesh_traversal cimport YTEmbreeScene</p>
<pre>from mesh_samplers cimport \
sample_hex, \
sample_tetra, \
sample_wedge</pre>
<p>-from pyembree.rtcore cimport \</p>
<ul><li><p>Vertex, \</p></li>
<li><p>Triangle, \</p></li>
<li><p>Vec3f</p></li></ul>
<pre>from mesh_intersection cimport \
patchIntersectFunc, \
patchBoundsFunc</pre>
<p>-from libc.stdlib cimport malloc, free -from libc.math cimport fmax, sqrt -cimport numpy as np</p>
<pre>from yt.utilities.exceptions import YTElementTypeNotRecognized
cdef extern from "mesh_construction.h":</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/3e0ff546abca/">https://bitbucket.org/yt_analysis/yt/commits/3e0ff546abca/</a> Changeset: 3e0ff546abca Branch: yt User: atmyers Date: 2016-04-24 03:29:08+00:00 Summary: some re-naming. Affected #: 2 files</p>
<p>diff -r 88ad574799ebf38c392490f4857d07763e9a7006 -r 3e0ff546abcaf0efb7347cb4db6a636844b130e7 yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -45,13 +45,13 @@</p>
<pre>cdef class BVH:
cdef BVHNode* root</pre>
<ul><li><p>cdef Triangle* triangles</p></li></ul>
<p>+ cdef void* primitives</p>
<pre>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</pre>
<ul><li><p>cdef np.int64_t num_tri</p></li></ul>
<p>+ cdef np.int64_t num_prim</p>
<pre>cdef np.int64_t num_elem
cdef np.int64_t num_verts_per_elem
cdef np.int64_t num_field_per_elem</pre>
<p>diff -r 88ad574799ebf38c392490f4857d07763e9a7006 -r 3e0ff546abcaf0efb7347cb4db6a636844b130e7 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -106,7 +106,7 @@</p>
<pre>self.num_tri_per_elem = TETRA_NT
tri_array = triangulate_tetra
self.sampler = P1Sampler</pre>
<ul><li><p>self.num_tri = self.num_tri_per_elem*self.num_elem</p></li></ul>
<p>+ self.num_prim = self.num_tri_per_elem*self.num_elem</p>
<pre># allocate storage
cdef np.int64_t v_size = self.num_verts_per_elem * self.num_elem * 3</pre>
<p>@@ -130,16 +130,16 @@</p>
<pre>cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri</pre>
<ul><li><p>self.triangles = <Triangle*> malloc(self.num_tri * sizeof(Triangle))</p></li>
<li><p>self.centroids = <np.float64_t**> malloc(self.num_tri * sizeof(np.float64_t*))</p></li>
<li><p>for i in range(self.num_tri):</p></li></ul>
<p>+ self.primitives = malloc(self.num_prim * sizeof(Triangle)) + self.centroids = <np.float64_t**> malloc(self.num_prim * sizeof(np.float64_t*)) + for i in range(self.num_prim):</p>
<pre>self.centroids[i] = <np.float64_t*> malloc(3*sizeof(np.float64_t))</pre>
<ul><li><p>self.bboxes = <BBox*> malloc(self.num_tri * sizeof(BBox))</p></li></ul>
<p>+ self.bboxes = <BBox*> malloc(self.num_prim * sizeof(BBox))</p>
<pre> 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</pre>
<ul><li><p>tri = &(self.triangles[tri_index])</p></li></ul>
<p>+ tri = &(<Triangle*> self.primitives)[tri_index]</p>
<pre>tri.elem_id = i
v0 = indices[i][tri_array[j][0]]
v1 = indices[i][tri_array[j][1]]</pre>
<p>@@ -148,14 +148,14 @@</p>
<pre>tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>self.get_centroid(self.triangles,</p></li></ul>
<p>+ self.get_centroid(self.primitives,</p>
<pre>tri_index,
self.centroids[tri_index])</pre>
<ul><li><p>self.get_bbox(self.triangles,</p></li></ul>
<p>+ self.get_bbox(self.primitives,</p>
<pre> tri_index,
&(self.bboxes[tri_index]))
</pre>
<ul><li><p>self.root = self._recursive_build(0, self.num_tri)</p></li></ul>
<p>+ self.root = self._recursive_build(0, self.num_prim)</p>
<pre>cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:</pre>
<p>@@ -165,9 +165,9 @@</p>
<pre>def __dealloc__(self):
self._recursive_free(self.root)</pre>
<ul><li><p>free(self.triangles)</p></li></ul>
<p>+ free(self.primitives)</p>
<pre>cdef np.int64_t i</pre>
<ul><li><p>for i in range(self.num_tri):</p></li></ul>
<p>+ for i in range(self.num_prim):</p>
<pre>free(self.centroids[i])
free(self.centroids)
free(self.bboxes)</pre>
<p>@@ -183,13 +183,15 @@</p>
<pre># 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
# will have centroids *greater* than "split" along "ax".</pre>
<p>+ cdef Triangle * triangles = <Triangle*> self.primitives +</p>
<pre> cdef np.int64_t mid = begin
while (begin != end):
if self.centroids[mid][ax] > split:
mid += 1
elif self.centroids[begin][ax] > split:</pre>
<ul><li><p>self.triangles[mid], self.triangles[begin] = \</p></li>
<li><p>self.triangles[begin], self.triangles[mid]</p></li></ul>
<p>+ triangles[mid], triangles[begin] = \ + triangles[begin], triangles[mid]</p>
<pre>self.centroids[mid], self.centroids[begin] = \
self.centroids[begin], self.centroids[mid]
self.bboxes[mid], self.bboxes[begin] = \</pre>
<p>@@ -255,7 +257,7 @@</p>
<pre> cdef Triangle* tri
if (node.end - node.begin) <= LEAF_SIZE:
for i in range(node.begin, node.end):</pre>
<ul><li><p>hit = self.get_intersect(self.triangles, i, ray)</p></li></ul>
<p>+ hit = self.get_intersect(self.primitives, i, ray)</p>
<pre> return
# if not leaf, intersect with left and right children</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/403971652a25/">https://bitbucket.org/yt_analysis/yt/commits/403971652a25/</a> Changeset: 403971652a25 Branch: yt User: atmyers Date: 2016-04-24 05:40:44+00:00 Summary: progress towards implementing second order for BVH. Affected #: 2 files</p>
<p>diff -r 3e0ff546abcaf0efb7347cb4db6a636844b130e7 -r 403971652a25420a13d9f61abd54aef0982e8ced yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -4,6 +4,20 @@</p>
<pre>from yt.utilities.lib.element_mappings cimport ElementSampler
from yt.utilities.lib.primitives cimport Triangle
</pre>
<p>+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] +</p>
<pre># ray data structure
cdef struct Ray:
np.float64_t origin[3]</pre>
<p>@@ -43,24 +57,33 @@</p>
<pre> const np.int64_t item,
BBox* bbox) nogil
</pre>
<p>+</p>
<pre>cdef class BVH:
cdef BVHNode* root
cdef void* primitives</pre>
<p>+ cdef np.int64_t* prim_ids</p>
<pre>cdef np.float64_t** centroids
cdef BBox* bboxes
cdef np.float64_t* vertices
cdef np.float64_t* field_data</pre>
<ul><li><p>cdef np.int64_t num_tri_per_elem</p></li></ul>
<p>+ cdef np.int64_t num_prim_per_elem</p>
<pre>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</pre>
<p>+ cdef int[MAX_NUM_TRI][3] tri_array</p>
<pre> 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</pre>
<p>+ 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</p>
<pre> cdef void intersect(self, Ray* ray) nogil
cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil</pre>
<p>diff -r 3e0ff546abcaf0efb7347cb4db6a636844b130e7 -r 403971652a25420a13d9f61abd54aef0982e8ced yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -9,35 +9,27 @@</p>
<pre>Triangle, \
ray_triangle_intersect, \
triangle_centroid, \</pre>
<ul><li><p>triangle_bbox</p></li></ul>
<p>+ triangle_bbox, \ + Patch, \ + ray_patch_intersect, \ + patch_centroid, \ + patch_bbox</p>
<pre>from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
P1Sampler3D, \</pre>
<ul><li><p>W1Sampler3D</p></li></ul>
<p>+ W1Sampler3D, \ + S2Sampler3D</p>
<pre>cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()</pre>
<p>+cdef ElementSampler S2Sampler = S2Sampler3D()</p>
<pre>cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
double fmin(double x, double y)
</pre>
<p>-cdef extern from “mesh_construction.h”:</p>
<ul><li><p>enum:</p></li>
<li><p>MAX_NUM_TRI</p></li></ul>
<p>–</p>
<ul><li><p>int HEX_NV</p></li>
<li><p>int HEX_NT</p></li>
<li><p>int TETRA_NV</p></li>
<li><p>int TETRA_NT</p></li>
<li><p>int WEDGE_NV</p></li>
<li><p>int WEDGE_NT</p></li>
<li><p>int triangulate_hex[MAX_NUM_TRI][3]</p></li>
<li><p>int triangulate_tetra[MAX_NUM_TRI][3]</p></li>
<li><p>int triangulate_wedge[MAX_NUM_TRI][3]</p></li></ul>
<p>–</p>
<pre># define some constants
cdef np.float64_t INF = np.inf
cdef np.int64_t LEAF_SIZE = 16</pre>
<p>@@ -88,34 +80,41 @@</p>
<pre> self.num_verts_per_elem = indices.shape[1]
self.num_field_per_elem = field_data.shape[1]
</pre>
<ul><li><p>self.get_centroid = triangle_centroid</p></li>
<li><p>self.get_bbox = triangle_bbox</p></li>
<li><p>self.get_intersect = ray_triangle_intersect</p></li></ul>
<p>–</p>
<pre># We need to figure out what kind of elements we've been handed.</pre>
<ul><li><p>cdef int[MAX_NUM_TRI][3] tri_array if self.num_verts_per_elem == 8:</p></li>
<li><p>self.num_tri_per_elem = HEX_NT</p></li>
<li><p>tri_array = triangulate_hex</p></li></ul>
<p>+ self.num_prim_per_elem = HEX_NT + self.tri_array = triangulate_hex</p>
<pre>self.sampler = Q1Sampler
elif self.num_verts_per_elem == 6:</pre>
<ul><li><p>self.num_tri_per_elem = WEDGE_NT</p></li>
<li><p>tri_array = triangulate_wedge</p></li></ul>
<p>+ self.num_prim_per_elem = WEDGE_NT + self.tri_array = triangulate_wedge</p>
<pre>self.sampler = W1Sampler
elif self.num_verts_per_elem == 4:</pre>
<ul><li><p>self.num_tri_per_elem = TETRA_NT</p></li>
<li><p>tri_array = triangulate_tetra</p></li></ul>
<p>+ self.num_prim_per_elem = TETRA_NT + self.tri_array = triangulate_tetra</p>
<pre>self.sampler = P1Sampler</pre>
<ul><li><p>self.num_prim = self.num_tri_per_elem*self.num_elem</p></li></ul>
<p>+ elif self.num_verts_per_elem == 20: + self.num_prim_per_elem = 8 + self.sampler = S2Sampler + else: + raise NotImplementedError + self.num_prim = self.num_prim_per_elem*self.num_elem</p>
<pre># 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))</pre>
<p>+ self.primitives = malloc(self.num_prim * sizeof(Triangle)) + 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))</p>
<pre># create data buffers</pre>
<ul><li><p>cdef np.int64_t i, j, k</p></li></ul>
<p>+ cdef np.int64_t j, k</p>
<pre> cdef np.int64_t field_offset, vertex_offset
for i in range(self.num_elem):
for j in range(self.num_verts_per_elem):</pre>
<p>@@ -126,24 +125,47 @@</p>
<pre> for j in range(self.num_field_per_elem):
self.field_data[field_offset + j] = field_data[i][j]
</pre>
<p>+ # set up primitives + if self.num_verts_per_elem == 20: + self.get_centroid = patch_centroid + self.get_bbox = patch_bbox + self.get_intersect = ray_patch_intersect + self._set_up_patches(vertices, indices) + else: + 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: + pass + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + cdef void _set_up_triangles(self, np.float64_t[:, :] vertices, + np.int64_t[:, :] indices) nogil:</p>
<pre># fill our array of primitives
cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri</pre>
<ul><li><p>self.primitives = malloc(self.num_prim * sizeof(Triangle))</p></li>
<li><p>self.centroids = <np.float64_t**> malloc(self.num_prim * sizeof(np.float64_t*))</p></li>
<li><p>for i in range(self.num_prim):</p></li>
<li><p>self.centroids[i] = <np.float64_t*> malloc(3*sizeof(np.float64_t))</p></li>
<li><p>self.bboxes = <BBox*> malloc(self.num_prim * sizeof(BBox))</p></li></ul>
<p>+ cdef np.int64_t i, j, k</p>
<pre>for i in range(self.num_elem):</pre>
<ul><li><p>offset = self.num_tri_per_elem*i</p></li>
<li><p>for j in range(self.num_tri_per_elem):</p></li></ul>
<p>+ offset = self.num_prim_per_elem*i + for j in range(self.num_prim_per_elem):</p>
<pre>tri_index = offset + j</pre>
<p>+ self.prim_ids[tri_index] = tri_index</p>
<pre>tri = &(<Triangle*> self.primitives)[tri_index]
tri.elem_id = i</pre>
<ul><li><p>v0 = indices[i][tri_array[j][0]]</p></li>
<li><p>v1 = indices[i][tri_array[j][1]]</p></li>
<li><p>v2 = indices[i][tri_array[j][2]]</p></li></ul>
<p>+ v0 = indices[i][self.tri_array[j][0]] + v1 = indices[i][self.tri_array[j][1]] + v2 = indices[i][self.tri_array[j][2]]</p>
<pre>for k in range(3):
tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]</pre>
<p>@@ -155,8 +177,6 @@</p>
<pre> tri_index,
&(self.bboxes[tri_index]))
</pre>
<ul><li><p>self.root = self._recursive_build(0, self.num_prim)</p></li></ul>
<p>–</p>
<pre> cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:
self._recursive_free(node.left)</pre>
<p>@@ -166,6 +186,7 @@</p>
<pre>def __dealloc__(self):
self._recursive_free(self.root)
free(self.primitives)</pre>
<p>+ free(self.prim_ids)</p>
<pre> cdef np.int64_t i
for i in range(self.num_prim):
free(self.centroids[i])</pre>
<p>@@ -179,19 +200,17 @@</p>
<pre> @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:</pre>
<ul><li><p># this re-orders the triangle array so that all of the triangles</p></li></ul>
<p>+ # this re-orders the primitive array so that all of the primitives</p>
<pre># to the left of mid have centroids less than or equal to "split"</pre>
<ul><li><p># along the direction “ax”. All the triangles to the right of mid</p></li></ul>
<p>+ # along the direction “ax”. All the primitives to the right of mid</p>
<pre># will have centroids *greater* than "split" along "ax".</pre>
<ul><li><p>cdef Triangle * triangles = <Triangle*> self.primitives</p></li></ul>
<p>–</p>
<pre> cdef np.int64_t mid = begin
while (begin != end):
if self.centroids[mid][ax] > split:
mid += 1
elif self.centroids[begin][ax] > split:</pre>
<ul><li><p>triangles[mid], triangles[begin] = \</p></li>
<li><p>triangles[begin], triangles[mid]</p></li></ul>
<p>+ self.prim_ids[mid], self.prim_ids[begin] = \ + self.prim_ids[begin], self.prim_ids[mid]</p>
<pre>self.centroids[mid], self.centroids[begin] = \
self.centroids[begin], self.centroids[mid]
self.bboxes[mid], self.bboxes[begin] = \</pre>
<p>@@ -254,10 +273,9 @@</p>
<pre># check for leaf
cdef np.int64_t i, hit</pre>
<ul><li><p>cdef Triangle* tri if (node.end – node.begin) <= LEAF_SIZE:</p>
<pre>for i in range(node.begin, node.end):</pre></li>
<li><p>hit = self.get_intersect(self.primitives, i, ray)</p></li></ul>
<p>+ hit = self.get_intersect(self.primitives, self.prim_ids[i], ray)</p>
<pre> return
# if not leaf, intersect with left and right children</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/65b621ce3d36/">https://bitbucket.org/yt_analysis/yt/commits/65b621ce3d36/</a> Changeset: 65b621ce3d36 Branch: yt User: atmyers Date: 2016-04-24 06:30:15+00:00 Summary: progress towards implementing second order for BVH. Affected #: 2 files</p>
<p>diff -r 403971652a25420a13d9f61abd54aef0982e8ced -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -17,6 +17,7 @@</p>
<pre>int triangulate_hex[MAX_NUM_TRI][3]
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]</pre>
<p>+ int hex20_faces[6][8]</p>
<pre># ray data structure
cdef struct Ray:</pre>
<p>diff -r 403971652a25420a13d9f61abd54aef0982e8ced -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -144,7 +144,28 @@</p>
<pre> @cython.cdivision(True)
cdef void _set_up_patches(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:</pre>
<ul><li><p>pass</p></li></ul>
<p>+ cdef Patch* patch + cdef np.ndarray[np.float64_t, ndim=2] element_vertices + 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 + element_vertices = vertices[indices[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 + 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] = element_vertices[ind][idim] + self.get_centroid(self.primitives, + prim_index, + self.centroids[prim_index]) + self.get_bbox(self.primitives, + prim_index, + &(self.bboxes[prim_index])) +</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)</pre>
<p>@@ -170,12 +191,12 @@</p>
<pre>tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>self.get_centroid(self.primitives,</p></li>
<li><p>tri_index,</p></li>
<li><p>self.centroids[tri_index])</p></li>
<li><p>self.get_bbox(self.primitives,</p></li>
<li><p>tri_index,</p></li>
<li><p>&(self.bboxes[tri_index]))</p></li></ul>
<p>+ self.get_centroid(self.primitives, + tri_index, + self.centroids[tri_index]) + self.get_bbox(self.primitives, + tri_index, + &(self.bboxes[tri_index]))</p>
<pre>cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/2f33070b9343/">https://bitbucket.org/yt_analysis/yt/commits/2f33070b9343/</a> Changeset: 2f33070b9343 Branch: yt User: atmyers Date: 2016-04-25 02:08:06+00:00 Summary: 20-node hexes working for BVH ray tracing. Affected #: 4 files</p>
<p>diff -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 -r 2f33070b93438f75f0992282dc96aefff8d5433f setup.py --- a/setup.py +++ b/setup.py @@ -171,6 +171,7 @@</p>
<pre>libraries=std_libs,
depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",
"yt/utilities/lib/element_mappings.pxd",</pre>
<p>+ “yt/utilities/lib/vec3_ops.pxd”,</p>
<pre> "yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],</pre>
<p>diff -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 -r 2f33070b93438f75f0992282dc96aefff8d5433f yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -20,6 +20,7 @@</p>
<pre>P1Sampler3D, \
W1Sampler3D, \
S2Sampler3D</pre>
<p>+from yt.utilities.lib.vec3_ops cimport L2_norm</p>
<pre>cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()</pre>
<p>@@ -94,7 +95,7 @@</p>
<pre>self.tri_array = triangulate_tetra
self.sampler = P1Sampler
elif self.num_verts_per_elem == 20:</pre>
<ul><li><p>self.num_prim_per_elem = 8</p></li></ul>
<p>+ self.num_prim_per_elem = 6</p>
<pre>self.sampler = S2Sampler
else:
raise NotImplementedError</pre>
<p>@@ -105,7 +106,6 @@</p>
<pre>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))</pre>
<ul><li><p>self.primitives = malloc(self.num_prim * sizeof(Triangle)) 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</p></li></ul>
<p>@@ -127,11 +127,13 @@</p>
<pre># set up primitives
if self.num_verts_per_elem == 20:</pre>
<p>+ self.primitives = malloc(self.num_prim * sizeof(Patch))</p>
<pre>self.get_centroid = patch_centroid
self.get_bbox = patch_bbox
self.get_intersect = ray_patch_intersect
self._set_up_patches(vertices, indices)
else:</pre>
<p>+ self.primitives = malloc(self.num_prim * sizeof(Triangle))</p>
<pre>self.get_centroid = triangle_centroid
self.get_bbox = triangle_bbox
self.get_intersect = ray_triangle_intersect</pre>
<p>@@ -145,27 +147,25 @@</p>
<pre> cdef void _set_up_patches(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:
cdef Patch* patch</pre>
<ul><li><p>cdef np.ndarray[np.float64_t, ndim=2] element_vertices cdef np.int64_t i, j, k, ind, idim cdef np.int64_t offset, prim_index for i in range(self.num_elem):</p>
<pre>offset = self.num_prim_per_elem*i</pre></li>
<li><p>element_vertices = vertices[indices[i]] for j in range(self.num_prim_per_elem): # for each face</p>
<pre>prim_index = offset + j
patch = &( <Patch*> self.primitives)[prim_index]
self.prim_ids[prim_index] = prim_index</pre></li></ul>
<p>+ patch.elem_id = i</p>
<pre>for k in range(8): # for each vertex
ind = hex20_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)</pre>
<ul><li><p>patch.v[k][idim] = element_vertices[ind][idim]</p></li></ul>
<p>+ patch.v[k][idim] = vertices[indices[i, ind]][idim]</p>
<pre>self.get_centroid(self.primitives,
prim_index,
self.centroids[prim_index])
self.get_bbox(self.primitives,
prim_index,
&(self.bboxes[prim_index]))</pre>
<p>–</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)</pre>
<p>@@ -265,9 +265,10 @@</p>
<pre> return
cdef np.float64_t[3] position</pre>
<p>+ cdef np.float64_t length = L2_norm(ray.direction)</p>
<pre>cdef np.int64_t i
for i in range(3):</pre>
<ul><li><p>position[i] = ray.origin[i] + ray.t_far*ray.direction[i]</p></li></ul>
<p>+ position[i] = ray.origin[i] + ray.t_far*ray.direction[i] / length</p>
<pre>cdef np.float64_t* vertex_ptr
cdef np.float64_t* field_ptr</pre>
<p>diff -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 -r 2f33070b93438f75f0992282dc96aefff8d5433f yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -4,7 +4,7 @@</p>
<pre>cimport cython.floating
from libc.math cimport fabs
</pre>
<p>-from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance +from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance, L2_norm</p>
<pre>from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox
cdef np.float64_t DETERMINANT_EPS = 1.0e-10</pre>
<p>diff -r 65b621ce3d369cd9c8bb0fb6c9be8ecc1500c291 -r 2f33070b93438f75f0992282dc96aefff8d5433f yt/visualization/volume_rendering/render_source.py --- a/yt/visualization/volume_rendering/render_source.py +++ b/yt/visualization/volume_rendering/render_source.py @@ -356,7 +356,7 @@</p>
<pre>self.field = field
self.volume = None
self.current_image = None</pre>
<ul><li><p>self.engine = ‘embree’</p></li></ul>
<p>+ self.engine = ‘bvh’</p>
<pre># default color map
self._cmap = ytcfg.get("yt", "default_colormap")</pre>
<p>@@ -491,15 +491,8 @@</p>
<pre> field_data = np.expand_dims(field_data, 1)
# Here, we decide whether to render based on high-order or</pre>
<ul><li><p># low-order geometry. Right now, high-order geometry is only</p></li>
<li><p># supported by Embree.</p></li>
<li><p>if indices.shape[1] == 20:</p></li>
<li><p># hexahedral</p></li>
<li><p>mylog.warning("20-node hexes not yet supported, " +</p></li>
<li><p>“dropping to 1st order.”)</p></li>
<li><p>field_data = field_data[:, 0:8]</p></li>
<li><p>indices = indices[:, 0:8]</p></li>
<li><p>elif indices.shape[1] == 27:</p></li></ul>
<p>+ # low-order geometry. + if indices.shape[1] == 27:</p>
<pre># hexahedral
mylog.warning("27-node hexes not yet supported, " +
"dropping to 1st order.")</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/c4fa0eadfddb/">https://bitbucket.org/yt_analysis/yt/commits/c4fa0eadfddb/</a> Changeset: c4fa0eadfddb Branch: yt User: atmyers Date: 2016-04-25 02:18:30+00:00 Summary: fixing tests. Affected #: 4 files</p>
<p>diff -r 2f33070b93438f75f0992282dc96aefff8d5433f -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -265,10 +265,9 @@</p>
<pre> return
cdef np.float64_t[3] position</pre>
<ul><li><p>cdef np.float64_t length = L2_norm(ray.direction) cdef np.int64_t i for i in range(3):</p></li>
<li><p>position[i] = ray.origin[i] + ray.t_far*ray.direction[i] / length</p></li></ul>
<p>+ position[i] = ray.origin[i] + ray.t_far*ray.direction[i]</p>
<pre>cdef np.float64_t* vertex_ptr
cdef np.float64_t* field_ptr</pre>
<p>diff -r 2f33070b93438f75f0992282dc96aefff8d5433f -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -209,7 +209,7 @@</p>
<pre> iterations += 1
# t is the distance along the ray to this hit</pre>
<ul><li><p>cdef np.float64_t t = distance(S, ray.origin)</p></li></ul>
<p>+ cdef np.float64_t t = distance(S, ray.origin) / L2_norm(ray.direction)</p>
<pre># only count this is it's the closest hit
if (t < ray.t_near or t > ray.t_far):</pre>
<p>diff -r 2f33070b93438f75f0992282dc96aefff8d5433f -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 yt/visualization/volume_rendering/render_source.py --- a/yt/visualization/volume_rendering/render_source.py +++ b/yt/visualization/volume_rendering/render_source.py @@ -356,7 +356,7 @@</p>
<pre>self.field = field
self.volume = None
self.current_image = None</pre>
<ul><li><p>self.engine = ‘bvh’</p></li></ul>
<p>+ self.engine = ‘embree’</p>
<pre># default color map
self._cmap = ytcfg.get("yt", "default_colormap")</pre>
<p>diff -r 2f33070b93438f75f0992282dc96aefff8d5433f -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 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 @@ -49,7 +49,7 @@</p>
<pre> return images
</pre>
<p>-def compare(ds, im, test_prefix, decimals=12): +def compare(ds, im, test_prefix, decimals=2):</p>
<pre> def mesh_render_image_func(filename_prefix):
return im.write_image(filename_prefix)
</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/3e2e9fb0afbc/">https://bitbucket.org/yt_analysis/yt/commits/3e2e9fb0afbc/</a> Changeset: 3e2e9fb0afbc Branch: yt User: atmyers Date: 2016-04-25 22:56:13+00:00 Summary: enable tests for both ray tracing engines. Affected #: 5 files</p>
<p>diff -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 tests/tests.yaml --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -44,7 +44,7 @@</p>
<pre>local_tipsy_000:
- yt/frontends/tipsy/tests/test_outputs.py
</pre>
<ul><li><p>local_varia_000:</p></li></ul>
<p>+ local_varia_001:</p>
<pre>- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py</pre>
<p>diff -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 yt/config.py --- a/yt/config.py +++ b/yt/config.py @@ -65,6 +65,7 @@</p>
<pre>chunk_size = '1000',
xray_data_dir = '/does/not/exist',
default_colormap = 'arbre',</pre>
<p>+ ray_tracing_engine = ‘embree’,</p>
<pre> )
# 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</pre>
<p>diff -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 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 @@</p>
<pre> from yt.utilities.lib import mesh_traversal
except ImportError:
mesh_traversal = NotAModule("pyembree")</pre>
<p>+ ytcfg["yt", “ray_tracing_engine”] = “yt”</p>
<pre>try:
from yt.utilities.lib import mesh_construction
except ImportError:
mesh_construction = NotAModule("pyembree")</pre>
<p>+ ytcfg["yt", “ray_tracing_engine”] = “yt”</p>
<pre>class RenderSource(ParallelAnalysisInterface):</pre>
<p>@@ -356,7 +358,7 @@</p>
<pre>self.field = field
self.volume = None
self.current_image = None</pre>
<ul><li><p>self.engine = ‘embree’</p></li></ul>
<p>+ self.engine = ytcfg.get("yt", “ray_tracing_engine”)</p>
<pre># default color map
self._cmap = ytcfg.get("yt", "default_colormap")</pre>
<p>@@ -374,11 +376,11 @@</p>
<pre> if self.engine == 'embree':
self.volume = mesh_traversal.YTEmbreeScene()
self.build_volume_embree()</pre>
<ul><li><p>elif self.engine == ‘bvh’:</p></li></ul>
<p>+ elif self.engine == ‘yt’:</p>
<pre>self.build_volume_bvh()
else:
raise NotImplementedError("Invalid ray-tracing engine selected. "</pre>
<ul><li><p>“Choices are ‘embree’ and 'bvh'.”)</p></li></ul>
<p>+ “Choices are ‘embree’ and 'yt'.”)</p>
<pre>def cmap():
'''</pre>
<p>diff -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 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,11 +23,19 @@</p>
<pre>MeshSource, \
Scene, \
create_scene</pre>
<p>+from yt.config import \ + ytcfg</p>
<p>+def compare(ds, im, test_prefix, decimals=12): + def mesh_render_image_func(filename_prefix): + return im.write_image(filename_prefix)</p>
<p>-@requires_module("pyembree") -def test_surface_mesh_render(): + mesh_render_image_func.__name__ = "func_{}".format(test_prefix) + test = GenericImageTest(ds, mesh_render_image_func, decimals) + test.prefix = test_prefix + return test</p>
<p>+def surface_mesh_render():</p>
<pre> images = []
ds = fake_tetrahedral_ds()</pre>
<p>@@ -47,75 +55,112 @@</p>
<pre> images.append(im)
return images</pre>
<p>+ +@requires_module("pyembree") +def test_surface_mesh_render_pyembree(): + ytcfg["yt", “ray_tracing_engine”] = “embree” + surface_mesh_render()</p>
<p>+def test_surface_mesh_render(): + ytcfg["yt", “ray_tracing_engine”] = “yt” + surface_mesh_render()</p>
<p>-def compare(ds, im, test_prefix, decimals=2):</p>
<ul><li><p>def mesh_render_image_func(filename_prefix):</p></li>
<li><p>return im.write_image(filename_prefix)</p></li></ul>
<p>–</p>
<ul><li><p>mesh_render_image_func.__name__ = "func_{}".format(test_prefix)</p></li>
<li><p>test = GenericImageTest(ds, mesh_render_image_func, decimals)</p></li>
<li><p>test.prefix = test_prefix</p></li>
<li><p>return test</p></li></ul>
<pre>hex8 = "MOOSE_sample_data/out.e-s010"
hex8_fields = [('connect1', 'diffused'), ('connect2', 'convected')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(hex8)
@requires_module("pyembree")</pre>
<p>+def test_hex8_render_pyembree(): + for field in hex8_fields: + yield hex8_render("embree", field) + +@requires_ds(hex8)</p>
<pre>def test_hex8_render():
for field in hex8_fields:</pre>
<ul><li><p>ds = data_dir_load(hex8, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_hex8_%s_%s” % field)</p></li></ul>
<p>+ yield hex8_render("yt", field)</p>
<pre>tet4 = "MOOSE_sample_data/high_order_elems_tet4_refine_out.e"
tet4_fields = [("connect1", "u")]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(tet4)
@requires_module("pyembree")</pre>
<p>+def test_tet4_render_pyembree(): + for field in tet4_fields: + yield tet4_render("embree", field) + +@requires_ds(tet4)</p>
<pre>def test_tet4_render():
for field in tet4_fields:</pre>
<ul><li><p>ds = data_dir_load(tet4, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_tet4_%s_%s” % field)</p></li></ul>
<p>+ yield tet4_render("yt", field)</p>
<pre>hex20 = "MOOSE_sample_data/mps_out.e"
hex20_fields = [('connect2', 'temp')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(hex20)
@requires_module("pyembree")</pre>
<p>+def test_hex20_render_pyembree(): + for field in hex20_fields: + yield hex20_render("embree", field) + +@requires_ds(hex20)</p>
<pre>def test_hex20_render():
for field in hex20_fields:</pre>
<ul><li><p>ds = data_dir_load(hex20, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_hex20_%s_%s” % field)</p></li></ul>
<p>+ yield hex20_render("yt", field)</p>
<pre>wedge6 = "MOOSE_sample_data/wedge_out.e"
wedge6_fields = [('connect1', 'diffused')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(wedge6)
@requires_module("pyembree")</pre>
<p>+def test_wedge6_render_pyembree(): + for field in wedge6_fields: + yield wedge6_render("embree", field) + +@requires_ds(wedge6)</p>
<pre>def test_wedge6_render():
for field in wedge6_fields:</pre>
<ul><li><p>ds = data_dir_load(wedge6, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_wedge6_%s_%s” % field)</p></li></ul>
<p>+ yield wedge6_render("yt", field)</p>
<p>– -@requires_ds(hex8) -@requires_module("pyembree") -def test_perspective_mesh_render(): +def perspective_mesh_render(engine): + ytcfg["yt", “ray_tracing_engine”] = engine</p>
<pre>ds = data_dir_load(hex8)
sc = create_scene(ds, ("connect2", "diffused"))</pre>
<p>–</p>
<pre>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')</pre>
<p>@@ -123,12 +168,19 @@</p>
<pre>cam.set_position(cam_pos, north_vector)
cam.resolution = (800, 800)
im = sc.render()</pre>
<ul><li><p>yield compare(ds, im, “perspective_mesh_render”)</p></li></ul>
<p>– + return compare(ds, im, “%s_perspective_mesh_render” % engine) + +@requires_ds(hex8) +@requires_module("pyembree") +def test_perspective_mesh_render_pyembree(): + yield perspective_mesh_render("embree")</p>
<pre>@requires_ds(hex8)</pre>
<p>-@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</p>
<pre>ds = data_dir_load(hex8)
sc = Scene()
cam = sc.add_camera(ds)</pre>
<p>@@ -137,12 +189,18 @@</p>
<pre>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)</pre>
<p>–</p>
<pre>ms1 = MeshSource(ds, ('connect1', 'diffused'))
ms2 = MeshSource(ds, ('connect2', 'diffused'))</pre>
<p>–</p>
<pre>sc.add_source(ms1)
sc.add_source(ms2)</pre>
<p>+ im = sc.render() + return compare(ds, im, “%s_composite_mesh_render” % engine) + +@requires_ds(hex8) +@requires_module("pyembree") +def test_composite_mesh_render_pyembree(): + yield composite_mesh_render("embree")</p>
<ul><li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “composite_mesh_render”)</p></li></ul>
<p>+@requires_ds(hex8) +def test_composite_mesh_render(): + yield composite_mesh_render("yt")</p>
<p>diff -r c4fa0eadfddb576a99b0534ba2cbba8ec0da1359 -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 yt/visualization/volume_rendering/utils.py --- a/yt/visualization/volume_rendering/utils.py +++ b/yt/visualization/volume_rendering/utils.py @@ -29,7 +29,7 @@</p>
<pre>kwargs = {'lens_type': params['lens_type']}
if engine == 'embree':
sampler = mesh_traversal.EmbreeMeshSampler(*args, **kwargs)</pre>
<ul><li><p>elif engine == ‘bvh’:</p></li></ul>
<p>+ elif engine == ‘yt’:</p>
<pre> sampler = mesh_traversal.BVHMeshSampler(*args, **kwargs)
return sampler
</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/3f125e11e1ab/">https://bitbucket.org/yt_analysis/yt/commits/3f125e11e1ab/</a> Changeset: 3f125e11e1ab Branch: yt User: atmyers Date: 2016-04-25 23:26:15+00:00 Summary: decrement version number to make sure tests are passing Affected #: 1 file</p>
<p>diff -r 3e2e9fb0afbc3db857ce488dfce51e0cc7ac6822 -r 3f125e11e1abd769be7fdf4d3f568b6aefad42a9 tests/tests.yaml --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -44,7 +44,7 @@</p>
<pre>local_tipsy_000:
- yt/frontends/tipsy/tests/test_outputs.py
</pre>
<ul><li><p>local_varia_001:</p></li></ul>
<p>+ local_varia_000:</p>
<pre>- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/e93b05586862/">https://bitbucket.org/yt_analysis/yt/commits/e93b05586862/</a> Changeset: e93b05586862 Branch: yt User: atmyers Date: 2016-04-26 04:43:20+00:00 Summary: fixing an error in fake_hexahedral_ds Affected #: 1 file</p>
<p>diff -r 3f125e11e1abd769be7fdf4d3f568b6aefad42a9 -r e93b05586862707336ff18f24023e0367c461b27 yt/testing.py --- a/yt/testing.py +++ b/yt/testing.py @@ -299,18 +299,16 @@</p>
<pre> from yt.frontends.stream.sample_data.hexahedral_mesh import \
_connectivity, _coordinates
</pre>
<ul><li><p>_connectivity -= 1 # this mesh has an offset of 1</p></li></ul>
<p>–</p>
<pre># the distance from the origin
node_data = {}
dist = np.sum(_coordinates**2, 1)</pre>
<ul><li><p>node_data[('connect1', 'test')] = dist[_connectivity]</p></li></ul>
<p>+ node_data[('connect1', 'test')] = dist[_connectivity-1]</p>
<pre> # each element gets a random number
elem_data = {}
elem_data[('connect1', 'elem')] = np.random.rand(_connectivity.shape[0])
</pre>
<ul><li><p>ds = load_unstructured_mesh(_connectivity,</p></li></ul>
<p>+ ds = load_unstructured_mesh(_connectivity-1,</p>
<pre>_coordinates,
node_data=node_data,
elem_data=elem_data)</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/46233ffd47c5/">https://bitbucket.org/yt_analysis/yt/commits/46233ffd47c5/</a> Changeset: 46233ffd47c5 Branch: yt User: atmyers Date: 2016-04-26 16:02:43+00:00 Summary: regenerate varia answers. Affected #: 1 file</p>
<p>diff -r e93b05586862707336ff18f24023e0367c461b27 -r 46233ffd47c5d5a5256f608f69872259334638e0 tests/tests.yaml --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -44,7 +44,7 @@</p>
<pre>local_tipsy_000:
- yt/frontends/tipsy/tests/test_outputs.py
</pre>
<ul><li><p>local_varia_000:</p></li></ul>
<p>+ local_varia_001:</p>
<pre>- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/a4d508a26437/">https://bitbucket.org/yt_analysis/yt/commits/a4d508a26437/</a> Changeset: a4d508a26437 Branch: yt User: atmyers Date: 2016-04-26 22:25:21+00:00 Summary: using fused types and a return type to avoid having two different patch intersect functions. Affected #: 3 files</p>
<p>diff -r 46233ffd47c5d5a5256f608f69872259334638e0 -r a4d508a26437ada288757f83a895c455d0a01aee yt/utilities/lib/mesh_intersection.pyx --- a/yt/utilities/lib/mesh_intersection.pyx +++ b/yt/utilities/lib/mesh_intersection.pyx @@ -21,10 +21,13 @@</p>
<pre>cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
from yt.utilities.lib.mesh_samplers cimport sample_hex20</pre>
<p>+from yt.utilities.lib.bounding_volume_hierarchy cimport BBox</p>
<pre>from yt.utilities.lib.primitives cimport \
patchSurfaceFunc, \
patchSurfaceDerivU, \</pre>
<ul><li><p>patchSurfaceDerivV</p></li></ul>
<p>+ patchSurfaceDerivV, \ + RayHitData, \ + compute_patch_hit</p>
<pre>from vec3_ops cimport dot, subtract, cross, distance
</pre>
<p>@@ -71,77 +74,20 @@</p>
<pre> cdef Patch patch = patches[item]
</pre>
<ul><li><p># first we compute the two planes that define the ray.</p></li>
<li><p>cdef float[3] n, N1, N2</p></li>
<li><p>cdef float A = dot(ray.dir, ray.dir)</p></li>
<li><p>for i in range(3):</p></li>
<li><p>n[i] = ray.dir[i] / A</p></li></ul>
<p>–</p>
<ul><li><p>if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):</p></li>
<li><p>N1[0] = n[1]</p></li>
<li><p>N1[1] =-n[0]</p></li>
<li><p>N1[2] = 0.0</p></li>
<li><p>else:</p></li>
<li><p>N1[0] = 0.0</p></li>
<li><p>N1[1] = n[2]</p></li>
<li><p>N1[2] =-n[1]</p></li>
<li><p>cross(N1, n, N2)</p></li></ul>
<p>–</p>
<ul><li><p>cdef float d1 = -dot(N1, ray.org)</p></li>
<li><p>cdef float d2 = -dot(N2, ray.org)</p></li></ul>
<p>–</p>
<ul><li><p># the initial guess is set to zero</p></li>
<li><p>cdef float u = 0.0</p></li>
<li><p>cdef float v = 0.0</p></li>
<li><p>cdef float[3] S</p></li>
<li><p>patchSurfaceFunc(patch.v, u, v, S)</p></li>
<li><p>cdef float fu = dot(N1, S) + d1</p></li>
<li><p>cdef float fv = dot(N2, S) + d2</p></li>
<li><p>cdef float err = fmax(fabs(fu), fabs(fv))</p></li></ul>
<p>–</p>
<ul><li><p># begin Newton interation</p></li>
<li><p>cdef float tol = 1.0e-5</p></li>
<li><p>cdef int iterations = 0</p></li>
<li><p>cdef int max_iter = 10</p></li>
<li><p>cdef float[3] Su</p></li>
<li><p>cdef float[3] Sv</p></li>
<li><p>cdef float J11, J12, J21, J22, det</p></li>
<li><p>while ((err > tol) and (iterations < max_iter)):</p></li>
<li><p># compute the Jacobian</p></li>
<li><p>patchSurfaceDerivU(patch.v, u, v, Su)</p></li>
<li><p>patchSurfaceDerivV(patch.v, u, v, Sv)</p></li>
<li><p>J11 = dot(N1, Su)</p></li>
<li><p>J12 = dot(N1, Sv)</p></li>
<li><p>J21 = dot(N2, Su)</p></li>
<li><p>J22 = dot(N2, Sv)</p></li>
<li><p>det = (J11*J22 – J12*J21)</p></li></ul>
<p>–</p>
<ul><li><p># update the u, v values</p></li>
<li><p>u -= ( J22*fu – J12*fv) / det</p></li>
<li><p>v -= (-J21*fu + J11*fv) / det</p></li></ul>
<p>–</p>
<ul><li><p>patchSurfaceFunc(patch.v, u, v, S)</p></li>
<li><p>fu = dot(N1, S) + d1</p></li>
<li><p>fv = dot(N2, S) + d2</p></li></ul>
<p>–</p>
<ul><li><p>err = fmax(fabs(fu), fabs(fv))</p></li>
<li><p>iterations += 1</p></li></ul>
<p>–</p>
<ul><li><p># t is the distance along the ray to this hit</p></li>
<li><p>cdef float t = distance(S, ray.org)</p></li></ul>
<p>+ cdef RayHitData hd = compute_patch_hit(patch.v, ray.org, ray.dir)</p>
<pre># only count this is it's the closest hit</pre>
<ul><li><p>if (t < ray.tnear or t > ray.Ng[0]):</p></li></ul>
<p>+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):</p>
<pre> return
</pre>
<ul><li><p>if (fabs(u) <= 1.0 and fabs(v) <= 1.0 and iterations < max_iter):</p></li></ul>
<p>+ if (fabs(hd.u) <= 1.0 and fabs(hd.v) <= 1.0 and hd.converged):</p>
<pre># we have a hit, so update ray information</pre>
<ul><li><p>ray.u = u</p></li>
<li><p>ray.v = v</p></li></ul>
<p>+ ray.u = hd.u + ray.v = hd.v</p>
<pre>ray.geomID = patch.geomID
ray.primID = item</pre>
<ul><li><p>ray.Ng[0] = t</p></li></ul>
<p>+ ray.Ng[0] = hd.t</p>
<pre># sample the solution at the calculated point
sample_hex20(patches, ray)</pre>
<p>diff -r 46233ffd47c5d5a5256f608f69872259334638e0 -r a4d508a26437ada288757f83a895c455d0a01aee yt/utilities/lib/primitives.pxd --- a/yt/utilities/lib/primitives.pxd +++ b/yt/utilities/lib/primitives.pxd @@ -1,9 +1,16 @@</p>
<pre>cimport cython</pre>
<p>+cimport cython.floating</p>
<pre>import numpy as np
cimport numpy as np
from vec3_ops cimport dot, subtract, cross
from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox
</pre>
<p>+cdef struct RayHitData: + np.float64_t u + np.float64_t v + np.float64_t t + np.int64_t converged +</p>
<pre>cdef struct Triangle:
np.float64_t p0[3]
np.float64_t p1[3]</pre>
<p>@@ -40,6 +47,10 @@</p>
<pre>const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil</pre>
<p>+ +cdef RayHitData compute_patch_hit(cython.floating[8][3] verts, + cython.floating[3] ray_origin, + cython.floating[3] ray_direction) nogil</p>
<pre>cdef inline np.int64_t ray_patch_intersect(const void* primitives,
const np.int64_t item,</pre>
<p>diff -r 46233ffd47c5d5a5256f608f69872259334638e0 -r a4d508a26437ada288757f83a895c455d0a01aee yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -1,4 +1,4 @@ -cimport cython +cimport cython</p>
<pre>import numpy as np
cimport numpy as np
cimport cython.floating</pre>
<p>@@ -13,6 +13,7 @@</p>
<pre> double fmax(double x, double y)
double fmin(double x, double y)
</pre>
<p>+</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>@@ -146,17 +147,15 @@</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef inline np.int64_t ray_patch_intersect(const void* primitives,</p>
<ul><li><p>const np.int64_t item,</p></li>
<li><p>Ray* ray) nogil:</p></li></ul>
<p>–</p>
<ul><li><p>cdef Patch patch = (<Patch*> primitives)[item]</p></li></ul>
<p>+cdef RayHitData compute_patch_hit(cython.floating[8][3] verts, + cython.floating[3] ray_origin, + cython.floating[3] ray_direction) nogil:</p>
<pre># first we compute the two planes that define the ray.</pre>
<ul><li><p>cdef np.float64_t[3] n, N1, N2</p></li>
<li><p>cdef np.float64_t A = dot(ray.direction, ray.direction)</p></li></ul>
<p>+ cdef cython.floating[3] n, N1, N2 + cdef cython.floating A = dot(ray_direction, ray_direction)</p>
<pre>for i in range(3):</pre>
<ul><li><p>n[i] = ray.direction[i] / A</p></li></ul>
<p>+ n[i] = ray_direction[i] / A</p>
<pre>if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):
N1[0] = n[1]</pre>
<p>@@ -168,29 +167,29 @@</p>
<pre> N1[2] =-n[1]
cross(N1, n, N2)
</pre>
<ul><li><p>cdef np.float64_t d1 = -dot(N1, ray.origin)</p></li>
<li><p>cdef np.float64_t d2 = -dot(N2, ray.origin)</p></li></ul>
<p>+ cdef cython.floating d1 = -dot(N1, ray_origin) + cdef cython.floating d2 = -dot(N2, ray_origin)</p>
<pre># the initial guess is set to zero</pre>
<ul><li><p>cdef np.float64_t u = 0.0</p></li>
<li><p>cdef np.float64_t v = 0.0</p></li>
<li><p>cdef np.float64_t[3] S</p></li>
<li><p>patchSurfaceFunc(patch.v, u, v, S)</p></li>
<li><p>cdef np.float64_t fu = dot(N1, S) + d1</p></li>
<li><p>cdef np.float64_t fv = dot(N2, S) + d2</p></li>
<li><p>cdef np.float64_t err = fmax(fabs(fu), fabs(fv))</p></li></ul>
<p>+ 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))</p>
<pre># begin Newton interation</pre>
<ul><li><p>cdef np.float64_t tol = 1.0e-5</p></li></ul>
<p>+ cdef cython.floating tol = 1.0e-5</p>
<pre>cdef int iterations = 0
cdef int max_iter = 10</pre>
<ul><li><p>cdef np.float64_t[3] Su</p></li>
<li><p>cdef np.float64_t[3] Sv</p></li>
<li><p>cdef np.float64_t J11, J12, J21, J22, det</p></li></ul>
<p>+ cdef cython.floating[3] Su + cdef cython.floating[3] Sv + cdef cython.floating J11, J12, J21, J22, det</p>
<pre>while ((err > tol) and (iterations < max_iter)):
# compute the Jacobian</pre>
<ul><li><p>patchSurfaceDerivU(patch.v, u, v, Su)</p></li>
<li><p>patchSurfaceDerivV(patch.v, u, v, Sv)</p></li></ul>
<p>+ patchSurfaceDerivU(verts, u, v, Su) + patchSurfaceDerivV(verts, u, v, Sv)</p>
<pre>J11 = dot(N1, Su)
J12 = dot(N1, Sv)
J21 = dot(N2, Su)</pre>
<p>@@ -201,7 +200,7 @@</p>
<pre>u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
</pre>
<ul><li><p>patchSurfaceFunc(patch.v, u, v, S)</p></li></ul>
<p>+ patchSurfaceFunc(verts, u, v, S)</p>
<pre> fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
</pre>
<p>@@ -209,16 +208,35 @@</p>
<pre> iterations += 1
# t is the distance along the ray to this hit</pre>
<ul><li><p>cdef np.float64_t t = distance(S, ray.origin) / L2_norm(ray.direction)</p></li></ul>
<p>+ 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</p>
<p>+ +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline np.int64_t ray_patch_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil: + + cdef Patch patch = (<Patch*> primitives)[item] + + cdef RayHitData hd = compute_patch_hit(patch.v, ray.origin, ray.direction) +</p>
<pre># only count this is it's the closest hit</pre>
<ul><li><p>if (t < ray.t_near or t > ray.t_far):</p></li></ul>
<p>+ if (hd.t < ray.t_near or hd.t > ray.t_far):</p>
<pre> return False
</pre>
<ul><li><p>if (fabs(u) <= 1.0 and fabs(v) <= 1.0 and iterations < max_iter):</p></li></ul>
<p>– + if (fabs(hd.u) <= 1.0 and fabs(hd.v) <= 1.0 and hd.converged):</p>
<pre># we have a hit, so update ray information</pre>
<ul><li><p>ray.t_far = t</p></li></ul>
<p>+ ray.t_far = hd.t</p>
<pre> ray.elem_id = patch.elem_id
return True
</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/371f59c79142/">https://bitbucket.org/yt_analysis/yt/commits/371f59c79142/</a> Changeset: 371f59c79142 Branch: yt User: atmyers Date: 2016-05-04 19:37:27+00:00 Summary: give a more informative error message when the element type can not be inferred. Affected #: 1 file</p>
<p>diff -r a4d508a26437ada288757f83a895c455d0a01aee -r 371f59c79142fe019351ffbd41bbe86e25b9d6cb yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -98,7 +98,8 @@</p>
<pre>self.num_prim_per_elem = 6
self.sampler = S2Sampler
else:</pre>
<ul><li><p>raise NotImplementedError</p></li></ul>
<p>+ raise NotImplementedError("Could not determine element type for " + "nverts = %d. " % self.num_verts_per_elem)</p>
<pre> self.num_prim = self.num_prim_per_elem*self.num_elem
# allocate storage</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/c3bd41ad7b51/">https://bitbucket.org/yt_analysis/yt/commits/c3bd41ad7b51/</a> Changeset: c3bd41ad7b51 Branch: yt User: atmyers Date: 2016-05-04 21:00:03+00:00 Summary: this type will be automatically inferred by Cython. Affected #: 1 file</p>
<p>diff -r 371f59c79142fe019351ffbd41bbe86e25b9d6cb -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -209,7 +209,6 @@</p>
<pre>self._recursive_free(self.root)
free(self.primitives)
free(self.prim_ids)</pre>
<ul><li><p>cdef np.int64_t i for i in range(self.num_prim):</p>
<pre>free(self.centroids[i])</pre>
<p>free(self.centroids)</p></li></ul>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/9f05c39bf3df/">https://bitbucket.org/yt_analysis/yt/commits/9f05c39bf3df/</a> Changeset: 9f05c39bf3df Branch: yt User: atmyers Date: 2016-05-04 21:24:56+00:00 Summary: moving Ray and BBox data structures and functions into primitives.pyx Affected #: 5 files</p>
<p>diff -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 -r 9f05c39bf3df7eae13e80815c667f8ba3c5cc9ed setup.py --- a/setup.py +++ b/setup.py @@ -169,8 +169,7 @@</p>
<pre>extra_compile_args=omp_args,
extra_link_args=omp_args,
libraries=std_libs,</pre>
<ul><li><p>depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",</p></li>
<li><p>“yt/utilities/lib/element_mappings.pxd”,</p></li></ul>
<p>+ depends=["yt/utilities/lib/element_mappings.pxd",</p>
<pre>"yt/utilities/lib/vec3_ops.pxd",
"yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",</pre>
<p>diff -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 -r 9f05c39bf3df7eae13e80815c667f8ba3c5cc9ed yt/utilities/lib/bounding_volume_hierarchy.pxd --- a/yt/utilities/lib/bounding_volume_hierarchy.pxd +++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd @@ -2,7 +2,7 @@</p>
<pre>import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler</pre>
<p>-from yt.utilities.lib.primitives cimport Triangle +from yt.utilities.lib.primitives cimport BBox, Ray</p>
<pre>cdef extern from "mesh_construction.h":
enum:</pre>
<p>@@ -19,22 +19,6 @@</p>
<pre> int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
</pre>
<p>-# ray data structure -cdef struct Ray:</p>
<ul><li><p>np.float64_t origin[3]</p></li>
<li><p>np.float64_t direction[3]</p></li>
<li><p>np.float64_t inv_dir[3]</p></li>
<li><p>np.float64_t data_val</p></li>
<li><p>np.float64_t t_near</p></li>
<li><p>np.float64_t t_far</p></li>
<li><p>np.int64_t elem_id</p></li>
<li><p>np.int64_t near_boundary</p></li></ul>
<p>– -# axis-aligned bounding box -cdef struct BBox:</p>
<ul><li><p>np.float64_t left_edge[3]</p></li>
<li><p>np.float64_t right_edge[3]</p></li></ul>
<p>–</p>
<pre># node for the bounding volume hierarchy
cdef struct BVHNode:
np.int64_t begin</pre>
<p>diff -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 -r 9f05c39bf3df7eae13e80815c667f8ba3c5cc9ed yt/utilities/lib/bounding_volume_hierarchy.pyx --- a/yt/utilities/lib/bounding_volume_hierarchy.pyx +++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx @@ -6,6 +6,9 @@</p>
<pre>from cython.parallel import parallel, prange
from yt.utilities.lib.primitives cimport \</pre>
<p>+ BBox, \ + Ray, \ + ray_bbox_intersect, \</p>
<pre>Triangle, \
ray_triangle_intersect, \
triangle_centroid, \</pre>
<p>@@ -36,26 +39,6 @@</p>
<pre>cdef np.int64_t LEAF_SIZE = 16
</pre>
<p>-@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil: -# <a href="https://tavianator.com/fast-branchless-raybounding-box-intersections/">https://tavianator.com/fast-branchless-raybounding-box-intersections/</a> –</p>
<ul><li><p>cdef np.float64_t tmin = -INF</p></li>
<li><p>cdef np.float64_t tmax = INF</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.int64_t i</p></li>
<li><p>cdef np.float64_t t1, t2</p></li>
<li><p>for i in range(3):</p></li>
<li><p>t1 = (bbox.left_edge[i] – ray.origin[i])*ray.inv_dir[i]</p></li>
<li><p>t2 = (bbox.right_edge[i] – ray.origin[i])*ray.inv_dir[i]</p></li>
<li><p>tmin = fmax(tmin, fmin(t1, t2))</p></li>
<li><p>tmax = fmin(tmax, fmax(t1, t2))</p></li></ul>
<p>–</p>
<ul><li><p>return tmax >= fmax(tmin, 0.0)</p></li></ul>
<p>– –</p>
<pre>cdef class BVH:
'''
</pre>
<p>diff -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 -r 9f05c39bf3df7eae13e80815c667f8ba3c5cc9ed yt/utilities/lib/primitives.pxd --- a/yt/utilities/lib/primitives.pxd +++ b/yt/utilities/lib/primitives.pxd @@ -3,7 +3,20 @@</p>
<pre>import numpy as np
cimport numpy as np
from vec3_ops cimport dot, subtract, cross</pre>
<p>-from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox + +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]</p>
<pre>cdef struct RayHitData:
np.float64_t u</pre>
<p>@@ -17,6 +30,8 @@</p>
<pre> np.float64_t p2[3]
np.int64_t elem_id
</pre>
<p>+cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil +</p>
<pre>cdef inline np.int64_t ray_triangle_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil</pre>
<p>diff -r c3bd41ad7b5151be1d053a31432ba74fdb3039c8 -r 9f05c39bf3df7eae13e80815c667f8ba3c5cc9ed yt/utilities/lib/primitives.pyx --- a/yt/utilities/lib/primitives.pyx +++ b/yt/utilities/lib/primitives.pyx @@ -5,14 +5,32 @@</p>
<pre>from libc.math cimport fabs
from yt.utilities.lib.vec3_ops cimport dot, subtract, cross, distance, L2_norm</pre>
<p>-from yt.utilities.lib.bounding_volume_hierarchy cimport Ray, BBox</p>
<pre>cdef np.float64_t DETERMINANT_EPS = 1.0e-10</pre>
<p>+cdef np.float64_t INF = np.inf</p>
<pre>cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
double fmin(double x, double y)
</pre>
<p>+@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil: +# <a href="https://tavianator.com/fast-branchless-raybounding-box-intersections/">https://tavianator.com/fast-branchless-raybounding-box-intersections/</a> + + 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)</p>
<pre>@cython.boundscheck(False)
@cython.wraparound(False)</pre>
<p><a href="https://bitbucket.org/yt_analysis/yt/commits/6068831a37f4/">https://bitbucket.org/yt_analysis/yt/commits/6068831a37f4/</a> Changeset: 6068831a37f4 Branch: yt User: xarthisius Date: 2016-05-11 20:04:07+00:00 Summary: Merged in atmyers/yt (pull request #2143)</p>
<p>High-Order Hexes for the Cython Ray Tracer Affected #: 15 files</p>
<p>diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 setup.py --- a/setup.py +++ b/setup.py @@ -157,14 +157,21 @@</p>
<pre> Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),</pre>
<p>+ 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"]),</p>
<pre> 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,</pre>
<ul><li><p>depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",</p></li>
<li><p>"yt/utilities/lib/vec3_ops.pxd"]),</p></li></ul>
<p>+ depends=["yt/utilities/lib/element_mappings.pxd", + “yt/utilities/lib/vec3_ops.pxd”, + "yt/utilities/lib/primitives.pxd"]),</p>
<pre> Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],
include_dirs=["yt/utilities/lib/",</pre>
<p>@@ -275,20 +282,30 @@</p>
<pre> embree_extensions = [
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],</pre>
<ul><li><p>depends=["yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ 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"]),</p>
<pre> Extension("yt.utilities.lib.mesh_traversal",
["yt/utilities/lib/mesh_traversal.pyx"],
depends=["yt/utilities/lib/mesh_traversal.pxd",</pre>
<ul><li><p>"yt/utilities/lib/grid_traversal.pxd"]),</p></li></ul>
<p>+ “yt/utilities/lib/grid_traversal.pxd”, + "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),</p>
<pre> 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",</pre>
<ul><li><p>"yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ “yt/utilities/lib/mesh_construction.pxd”, + “yt/utilities/lib/bounding_volume_hierarchy.pxd”, + "yt/utilities/lib/primitives.pxd"]),</p>
<pre> Extension("yt.utilities.lib.mesh_intersection",
["yt/utilities/lib/mesh_intersection.pyx"],
depends=["yt/utilities/lib/mesh_intersection.pxd",</pre>
<ul><li><p>"yt/utilities/lib/mesh_construction.pxd"]),</p></li></ul>
<p>+ “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"]),</p>
<pre> ]
embree_prefix = os.path.abspath(read_embree_location())</pre>
<p>diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 tests/tests.yaml --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -47,7 +47,7 @@</p>
<pre>local_tipsy_000:
- yt/frontends/tipsy/tests/test_outputs.py
</pre>
<ul><li><p>local_varia_000:</p></li></ul>
<p>+ local_varia_001:</p>
<pre>- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py</pre>
<p>diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/config.py --- a/yt/config.py +++ b/yt/config.py @@ -65,6 +65,7 @@</p>
<pre>chunk_size = '1000',
xray_data_dir = '/does/not/exist',
default_colormap = 'arbre',</pre>
<p>+ ray_tracing_engine = ‘embree’,</p>
<pre> )
# 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</pre>
<p>diff -r 5ffd0fbe9057407dbb37fa2225906315695d569a -r 6068831a37f4d92db1f3788c867e56efb286e587 yt/testing.py --- a/yt/testing.py +++ b/yt/testing.py @@ -306,18 +306,16 @@</p>
<pre> from yt.frontends.stream.sample_data.hexahedral_mesh import \
_connectivity, _coordinates
</pre>
<ul><li><p>_connectivity -= 1 # this mesh has an offset of 1</p></li></ul>
<p>–</p>
<pre># the distance from the origin
node_data = {}
dist = np.sum(_coordinates**2, 1)</pre>
<ul><li><p>node_data[('connect1', 'test')] = dist[_connectivity]</p></li></ul>
<p>+ node_data[('connect1', 'test')] = dist[_connectivity-1]</p>
<pre> # each element gets a random number
elem_data = {}
elem_data[('connect1', 'elem')] = np.random.rand(_connectivity.shape[0])
</pre>
<ul><li><p>ds = load_unstructured_mesh(_connectivity,</p></li></ul>
<p>+ ds = load_unstructured_mesh(_connectivity-1,</p>
<pre>_coordinates,
node_data=node_data,
elem_data=elem_data)</pre>
<p>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 @@</p>
<pre>import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler</pre>
<p>+from yt.utilities.lib.primitives cimport BBox, Ray</p>
<p>-# ray data structure -cdef struct Ray:</p>
<ul><li><p>np.float64_t origin[3]</p></li>
<li><p>np.float64_t direction[3]</p></li>
<li><p>np.float64_t inv_dir[3]</p></li>
<li><p>np.float64_t data_val</p></li>
<li><p>np.float64_t t_near</p></li>
<li><p>np.float64_t t_far</p></li>
<li><p>np.int64_t elem_id</p></li>
<li><p>np.int64_t near_boundary</p></li></ul>
<p>+cdef extern from “mesh_construction.h”: + enum: + MAX_NUM_TRI</p>
<p>-# axis-aligned bounding box -cdef struct BBox:</p>
<ul><li><p>np.float64_t left_edge[3]</p></li>
<li><p>np.float64_t right_edge[3]</p></li></ul>
<p>+ 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]</p>
<pre># node for the bounding volume hierarchy
cdef struct BVHNode:</pre>
<p>@@ -26,29 +26,49 @@</p>
<pre>BVHNode* left
BVHNode* right
BBox bbox</pre>
<p>– -# triangle data structure -cdef struct Triangle:</p>
<ul><li><p>np.float64_t p0[3]</p></li>
<li><p>np.float64_t p1[3]</p></li>
<li><p>np.float64_t p2[3]</p></li>
<li><p>np.int64_t elem_id</p></li>
<li><p>np.float64_t centroid[3]</p></li>
<li><p>BBox bbox</p></li></ul>
<p>+ +# 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 <strong>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</strong> bbox) nogil +</p>
<pre>cdef class BVH:
cdef BVHNode* root</pre>
<ul><li><p>cdef Triangle* triangles</p></li></ul>
<p>+ cdef void* primitives + cdef np.int64_t* prim_ids + cdef np.float64_t** centroids + cdef BBox* bboxes</p>
<pre>cdef np.float64_t* vertices
cdef np.float64_t* field_data</pre>
<ul><li><p>cdef np.int64_t num_tri_per_elem</p></li>
<li><p>cdef np.int64_t num_tri</p></li></ul>
<p>+ cdef np.int64_t num_prim_per_elem + cdef np.int64_t num_prim</p>
<pre>cdef np.int64_t num_elem
cdef np.int64_t num_verts_per_elem
cdef np.int64_t num_field_per_elem</pre>
<p>+ cdef int[MAX_NUM_TRI][3] tri_array</p>
<pre>cdef ElementSampler sampler</pre>
<p>+ cdef centroid_func_type get_centroid + cdef bbox_func_type get_bbox + cdef intersect_func_type get_intersect</p>
<pre> cdef np.int64_t _partition(self, np.int64_t begin, np.int64_t end,
np.int64_t ax, np.float64_t split) nogil</pre>
<p>+ 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</p>
<pre> cdef void intersect(self, Ray* ray) nogil
cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil</pre>
<p>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 @@</p>
<pre>from libc.math cimport fabs
from libc.stdlib cimport malloc, free
from cython.parallel import parallel, prange</pre>
<p>-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</p>
<pre>from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
P1Sampler3D, \</pre>
<ul><li><p>W1Sampler3D</p></li></ul>
<p>+ W1Sampler3D, \ + S2Sampler3D +from yt.utilities.lib.vec3_ops cimport L2_norm</p>
<pre>cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()</pre>
<p>+cdef ElementSampler S2Sampler = S2Sampler3D()</p>
<pre>cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
double fmin(double x, double y)
</pre>
<p>-cdef extern from “mesh_construction.h”:</p>
<ul><li><p>enum:</p></li>
<li><p>MAX_NUM_TRI</p></li></ul>
<p>–</p>
<ul><li><p>int HEX_NV</p></li>
<li><p>int HEX_NT</p></li>
<li><p>int TETRA_NV</p></li>
<li><p>int TETRA_NT</p></li>
<li><p>int WEDGE_NV</p></li>
<li><p>int WEDGE_NT</p></li>
<li><p>int triangulate_hex[MAX_NUM_TRI][3]</p></li>
<li><p>int triangulate_tetra[MAX_NUM_TRI][3]</p></li>
<li><p>int triangulate_wedge[MAX_NUM_TRI][3]</p></li></ul>
<p>–</p>
<pre># define some constants</pre>
<p>-cdef np.float64_t DETERMINANT_EPS = 1.0e-10</p>
<pre>cdef np.float64_t INF = np.inf
cdef np.int64_t LEAF_SIZE = 16
</pre>
<p>-@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef np.int64_t ray_triangle_intersect(Ray* ray, const Triangle* tri) nogil: -# <a href="https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm">https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm</a> –</p>
<ul><li><p># edge vectors</p></li>
<li><p>cdef np.float64_t e1[3]</p></li>
<li><p>cdef np.float64_t e2[3]</p></li>
<li><p>subtract(tri.p1, tri.p0, e1)</p></li>
<li><p>subtract(tri.p2, tri.p0, e2)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t P[3]</p></li>
<li><p>cross(ray.direction, e2, P)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t det, inv_det</p></li>
<li><p>det = dot(e1, P)</p></li>
<li><p>if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):</p></li>
<li><p>return False</p></li>
<li><p>inv_det = 1.0 / det</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t T[3]</p></li>
<li><p>subtract(ray.origin, tri.p0, T)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t u = dot(T, P) * inv_det</p></li>
<li><p>if(u < 0.0 or u > 1.0):</p></li>
<li><p>return False</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t Q[3]</p></li>
<li><p>cross(T, e1, Q)</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t v = dot(ray.direction, Q) * inv_det</p></li>
<li><p>if(v < 0.0 or u + v > 1.0):</p></li>
<li><p>return False</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.float64_t t = dot(e2, Q) * inv_det</p></li></ul>
<p>–</p>
<ul><li><p>if(t > DETERMINANT_EPS and t < ray.t_far):</p></li>
<li><p>ray.t_far = t</p></li>
<li><p>ray.elem_id = tri.elem_id</p></li>
<li><p>return True</p></li></ul>
<p>–</p>
<ul><li><p>return False</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil: -# <a href="https://tavianator.com/fast-branchless-raybounding-box-intersections/">https://tavianator.com/fast-branchless-raybounding-box-intersections/</a> –</p>
<ul><li><p>cdef np.float64_t tmin = -INF</p></li>
<li><p>cdef np.float64_t tmax = INF</p></li></ul>
<p>–</p>
<ul><li><p>cdef np.int64_t i</p></li>
<li><p>cdef np.float64_t t1, t2</p></li>
<li><p>for i in range(3):</p></li>
<li><p>t1 = (bbox.left_edge[i] – ray.origin[i])*ray.inv_dir[i]</p></li>
<li><p>t2 = (bbox.right_edge[i] – ray.origin[i])*ray.inv_dir[i]</p></li>
<li><p>tmin = fmax(tmin, fmin(t1, t2))</p></li>
<li><p>tmax = fmin(tmax, fmax(t1, t2))</p></li></ul>
<p>–</p>
<ul><li><p>return tmax >= fmax(tmin, 0.0)</p></li></ul>
<p>– –</p>
<pre>cdef class BVH:
'''
</pre>
<p>@@ -130,29 +65,40 @@</p>
<pre> self.num_field_per_elem = field_data.shape[1]
# We need to figure out what kind of elements we've been handed.</pre>
<ul><li><p>cdef int[MAX_NUM_TRI][3] tri_array if self.num_verts_per_elem == 8:</p></li>
<li><p>self.num_tri_per_elem = HEX_NT</p></li>
<li><p>tri_array = triangulate_hex</p></li></ul>
<p>+ self.num_prim_per_elem = HEX_NT + self.tri_array = triangulate_hex</p>
<pre>self.sampler = Q1Sampler
elif self.num_verts_per_elem == 6:</pre>
<ul><li><p>self.num_tri_per_elem = WEDGE_NT</p></li>
<li><p>tri_array = triangulate_wedge</p></li></ul>
<p>+ self.num_prim_per_elem = WEDGE_NT + self.tri_array = triangulate_wedge</p>
<pre>self.sampler = W1Sampler
elif self.num_verts_per_elem == 4:</pre>
<ul><li><p>self.num_tri_per_elem = TETRA_NT</p></li>
<li><p>tri_array = triangulate_tetra</p></li></ul>
<p>+ self.num_prim_per_elem = TETRA_NT + self.tri_array = triangulate_tetra</p>
<pre>self.sampler = P1Sampler</pre>
<ul><li><p>self.num_tri = self.num_tri_per_elem*self.num_elem</p></li></ul>
<p>+ 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</p>
<pre># 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))</pre>
<p>+ 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))</p>
<pre># create data buffers</pre>
<ul><li><p>cdef np.int64_t i, j, k</p></li></ul>
<p>+ cdef np.int64_t j, k</p>
<pre> cdef np.int64_t field_offset, vertex_offset
for i in range(self.num_elem):
for j in range(self.num_verts_per_elem):</pre>
<p>@@ -163,29 +109,78 @@</p>
<pre> for j in range(self.num_field_per_elem):
self.field_data[field_offset + j] = field_data[i][j]
</pre>
<ul><li><p># fill our array of triangles</p></li></ul>
<p>+ # 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</p>
<pre>cdef np.int64_t offset, tri_index
cdef np.int64_t v0, v1, v2
cdef Triangle* tri</pre>
<ul><li><p>self.triangles = <Triangle*> malloc(self.num_tri * sizeof(Triangle))</p></li></ul>
<p>+ cdef np.int64_t i, j, k</p>
<pre>for i in range(self.num_elem):</pre>
<ul><li><p>offset = self.num_tri_per_elem*i</p></li>
<li><p>for j in range(self.num_tri_per_elem):</p></li></ul>
<p>+ offset = self.num_prim_per_elem*i + for j in range(self.num_prim_per_elem):</p>
<pre>tri_index = offset + j</pre>
<ul><li><p>tri = &(self.triangles[tri_index])</p></li></ul>
<p>+ self.prim_ids[tri_index] = tri_index + tri = &(<Triangle*> self.primitives)[tri_index]</p>
<pre>tri.elem_id = i</pre>
<ul><li><p>v0 = indices[i][tri_array[j][0]]</p></li>
<li><p>v1 = indices[i][tri_array[j][1]]</p></li>
<li><p>v2 = indices[i][tri_array[j][2]]</p></li></ul>
<p>+ v0 = indices[i][self.tri_array[j][0]] + v1 = indices[i][self.tri_array[j][1]] + v2 = indices[i][self.tri_array[j][2]]</p>
<pre>for k in range(3):
tri.p0[k] = vertices[v0][k]
tri.p1[k] = vertices[v1][k]
tri.p2[k] = vertices[v2][k]</pre>
<ul><li><p>tri.centroid[k] = (tri.p0[k] + tri.p1[k] + tri.p2[k]) / 3.0</p></li>
<li><p>tri.bbox.left_edge[k] = fmin(fmin(tri.p0[k], tri.p1[k]), tri.p2[k])</p></li>
<li><p>tri.bbox.right_edge[k] = fmax(fmax(tri.p0[k], tri.p1[k]), tri.p2[k])</p></li></ul>
<p>–</p>
<ul><li><p>self.root = self._recursive_build(0, self.num_tri)</p></li></ul>
<p>+ self.get_centroid(self.primitives, + tri_index, + self.centroids[tri_index]) + self.get_bbox(self.primitives, + tri_index, + &(self.bboxes[tri_index]))</p>
<pre>cdef void _recursive_free(self, BVHNode* node) nogil:
if node.end - node.begin > LEAF_SIZE:</pre>
<p>@@ -195,7 +190,12 @@</p>
<pre>def __dealloc__(self):
self._recursive_free(self.root)</pre>
<ul><li><p>free(self.triangles)</p></li></ul>
<p>+ free(self.primitives) + free(self.prim_ids) + for i in range(self.num_prim): + free(self.centroids[i]) + free(self.centroids) + free(self.bboxes)</p>
<pre> free(self.field_data)
free(self.vertices)
</pre>
<p>@@ -204,17 +204,21 @@</p>
<pre> @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:</pre>
<ul><li><p># this re-orders the triangle array so that all of the triangles</p></li>
<li><p># to the left of mid have centroids less than or equal to “split”</p></li>
<li><p># along the direction “ax”. All the triangles to the right of mid</p></li></ul>
<p>+ # 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</p>
<pre># will have centroids *greater* than "split" along "ax".
cdef np.int64_t mid = begin
while (begin != end):</pre>
<ul><li><p>if self.triangles[mid].centroid[ax] > split:</p></li></ul>
<p>+ if self.centroids[mid][ax] > split:</p>
<pre>mid += 1</pre>
<ul><li><p>elif self.triangles[begin].centroid[ax] > split:</p></li>
<li><p>self.triangles[mid], self.triangles[begin] = \</p></li>
<li><p>self.triangles[begin], self.triangles[mid]</p></li></ul>
<p>+ 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]</p>
<pre> mid += 1
begin += 1
return mid</pre>
<p>@@ -225,13 +229,13 @@</p>
<pre> cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil:
cdef np.int64_t i, j</pre>
<ul><li><p>cdef BBox box = self.triangles[begin].bbox</p></li></ul>
<p>+ cdef BBox box = self.bboxes[begin]</p>
<pre> for i in range(begin+1, end):
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],</pre>
<ul><li><p>self.triangles[i].bbox.left_edge[j])</p></li></ul>
<p>+ self.bboxes[i].left_edge[j])</p>
<pre>box.right_edge[j] = fmax(box.right_edge[j],</pre>
<ul><li><p>self.triangles[i].bbox.right_edge[j])</p></li></ul>
<p>+ self.bboxes[i].right_edge[j])</p>
<pre> node.bbox = box
@cython.boundscheck(False)</pre>
<p>@@ -249,7 +253,7 @@</p>
<pre>position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
cdef np.float64_t* vertex_ptr</pre>
<ul><li><p>cdef np.float64_t* field_ptr</p></li></ul>
<p>+ cdef np.float64_t* field_ptr</p>
<pre> 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
</pre>
<p>@@ -273,11 +277,9 @@</p>
<pre># check for leaf
cdef np.int64_t i, hit</pre>
<ul><li><p>cdef Triangle* tri if (node.end – node.begin) <= LEAF_SIZE:</p>
<pre>for i in range(node.begin, node.end):</pre></li>
<li><p>tri = &(self.triangles[i])</p></li>
<li><p>hit = ray_triangle_intersect(ray, tri)</p></li></ul>
<p>+ hit = self.get_intersect(self.primitives, self.prim_ids[i], ray)</p>
<pre> return
# if not leaf, intersect with left and right children</pre>
<p>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 @@</p>
<pre>import numpy as np
cimport cython</pre>
<p>+from libc.stdlib cimport malloc, free +from libc.math cimport fmax, sqrt +cimport numpy as np +</p>
<pre>cimport pyembree.rtcore as rtc</pre>
<p>-from mesh_traversal cimport YTEmbreeScene</p>
<pre>cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu</pre>
<p>+from pyembree.rtcore cimport \ + Vertex, \ + Triangle, \ + Vec3f + +from mesh_traversal cimport YTEmbreeScene</p>
<pre>from mesh_samplers cimport \
sample_hex, \
sample_tetra, \
sample_wedge</pre>
<p>-from pyembree.rtcore cimport \</p>
<ul><li><p>Vertex, \</p></li>
<li><p>Triangle, \</p></li>
<li><p>Vec3f</p></li></ul>
<pre>from mesh_intersection cimport \
patchIntersectFunc, \
patchBoundsFunc</pre>
<p>-from libc.stdlib cimport malloc, free -from libc.math cimport fmax, sqrt -cimport numpy as np</p>
<pre>from yt.utilities.exceptions import YTElementTypeNotRecognized
cdef extern from "mesh_construction.h":</pre>
<p>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 @@</p>
<pre>from yt.utilities.lib.mesh_construction cimport Patch
cimport cython
</pre>
<p>-cdef void patchIntersectFunc(Patch* patches,</p>
<ul><li><p>rtcr.RTCRay& ray,</p></li></ul>
<p>+cdef void patchIntersectFunc(Patch* patches, + rtcr.RTCRay& ray,</p>
<pre> size_t item) nogil
</pre>
<p>-cdef void patchBoundsFunc(Patch* patches,</p>
<ul><li><p>size_t item,</p></li></ul>
<p>+cdef void patchBoundsFunc(Patch* patches, + size_t item,</p>
<pre>rtcg.RTCBounds* bounds_o) nogil</pre>
<p>– -cdef void patchSurfaceFunc(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] S) nogil</p></li></ul>
<p>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 @@</p>
<pre>cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
from yt.utilities.lib.mesh_samplers cimport sample_hex20</pre>
<p>+from yt.utilities.lib.bounding_volume_hierarchy cimport BBox +from yt.utilities.lib.primitives cimport \ + patchSurfaceFunc, \ + patchSurfaceDerivU, \ + patchSurfaceDerivV, \ + RayHitData, \ + compute_patch_hit</p>
<pre>from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)</pre>
<p>-cdef void patchSurfaceFunc(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] S) nogil:</p></li></ul>
<p>–</p>
<ul><li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>S[i] = 0.25*(1.0 – u)*(1.0 – v)*(-u – v – 1)*patch.v[0][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 – v)*( u – v – 1)*patch.v[1][i] + \</p></li>
<li><p>0.25*(1.0 + u)*(1.0 + v)*( u + v – 1)*patch.v[2][i] + \</p></li>
<li><p>0.25*(1.0 – u)*(1.0 + v)*(-u + v – 1)*patch.v[3][i] + \</p></li>
<li><p>0.5*(1 – u)*(1 – v*v)*patch.v[4][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 – v)*patch.v[5][i] + \</p></li>
<li><p>0.5*(1 + u)*(1 – v*v)*patch.v[6][i] + \</p></li>
<li><p>0.5*(1 – u*u)*(1 + v)*patch.v[7][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void patchSurfaceDerivU(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] Su) nogil:</p></li>
<li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>Su[i] = (-0.25*(v – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(v – 1.0)*(u – v – 1) – 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(v + 1.0)*(u – v + 1) + 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[4][i] + u*(v – 1.0)*patch.v[5][i] – \</p></li>
<li><p>0.5*(v*v – 1.0)*patch.v[6][i] – u*(v + 1.0)*patch.v[7][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True) -cdef void patchSurfaceDerivV(const Patch& patch,</p>
<ul><li><p>const float u,</p></li>
<li><p>const float v,</p></li>
<li><p>float[3] Sv) nogil:</p></li>
<li><p>cdef int i</p></li>
<li><p>for i in range(3):</p></li>
<li><p>Sv[i] = (-0.25*(u – 1.0)*(u + v + 1) – 0.25*(u – 1.0)*(v – 1.0))*patch.v[0][i] + \</p></li>
<li><p>(-0.25*(u + 1.0)*(u – v – 1) + 0.25*(u + 1.0)*(v – 1.0))*patch.v[1][i] + \</p></li>
<li><p>( 0.25*(u + 1.0)*(u + v – 1) + 0.25*(u + 1.0)*(v + 1.0))*patch.v[2][i] + \</p></li>
<li><p>( 0.25*(u – 1.0)*(u – v + 1) – 0.25*(u – 1.0)*(v + 1.0))*patch.v[3][i] + \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[5][i] + v*(u – 1.0)*patch.v[4][i] – \</p></li>
<li><p>0.5*(u*u – 1.0)*patch.v[7][i] – v*(u + 1.0)*patch.v[6][i]</p></li></ul>
<p>– – -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.cdivision(True)</p>
<pre>cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:</pre>
<p>@@ -121,77 +74,20 @@</p>
<pre> cdef Patch patch = patches[item]
</pre>
<ul><li><p># first we compute the two planes that define the ray.</p></li>
<li><p>cdef float[3] n, N1, N2</p></li>
<li><p>cdef float A = dot(ray.dir, ray.dir)</p></li>
<li><p>for i in range(3):</p></li>
<li><p>n[i] = ray.dir[i] / A</p></li></ul>
<p>–</p>
<ul><li><p>if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):</p></li>
<li><p>N1[0] = n[1]</p></li>
<li><p>N1[1] =-n[0]</p></li>
<li><p>N1[2] = 0.0</p></li>
<li><p>else:</p></li>
<li><p>N1[0] = 0.0</p></li>
<li><p>N1[1] = n[2]</p></li>
<li><p>N1[2] =-n[1]</p></li>
<li><p>cross(N1, n, N2)</p></li></ul>
<p>–</p>
<ul><li><p>cdef float d1 = -dot(N1, ray.org)</p></li>
<li><p>cdef float d2 = -dot(N2, ray.org)</p></li></ul>
<p>–</p>
<ul><li><p># the initial guess is set to zero</p></li>
<li><p>cdef float u = 0.0</p></li>
<li><p>cdef float v = 0.0</p></li>
<li><p>cdef float[3] S</p></li>
<li><p>patchSurfaceFunc(patch, u, v, S)</p></li>
<li><p>cdef float fu = dot(N1, S) + d1</p></li>
<li><p>cdef float fv = dot(N2, S) + d2</p></li>
<li><p>cdef float err = fmax(fabs(fu), fabs(fv))</p></li></ul>
<p>–</p>
<ul><li><p># begin Newton interation</p></li>
<li><p>cdef float tol = 1.0e-5</p></li>
<li><p>cdef int iterations = 0</p></li>
<li><p>cdef int max_iter = 10</p></li>
<li><p>cdef float[3] Su</p></li>
<li><p>cdef float[3] Sv</p></li>
<li><p>cdef float J11, J12, J21, J22, det</p></li>
<li><p>while ((err > tol) and (iterations < max_iter)):</p></li>
<li><p># compute the Jacobian</p></li>
<li><p>patchSurfaceDerivU(patch, u, v, Su)</p></li>
<li><p>patchSurfaceDerivV(patch, u, v, Sv)</p></li>
<li><p>J11 = dot(N1, Su)</p></li>
<li><p>J12 = dot(N1, Sv)</p></li>
<li><p>J21 = dot(N2, Su)</p></li>
<li><p>J22 = dot(N2, Sv)</p></li>
<li><p>det = (J11*J22 – J12*J21)</p></li></ul>
<p>–</p>
<ul><li><p># update the u, v values</p></li>
<li><p>u -= ( J22*fu – J12*fv) / det</p></li>
<li><p>v -= (-J21*fu + J11*fv) / det</p></li></ul>
<p>–</p>
<ul><li><p>patchSurfaceFunc(patch, u, v, S)</p></li>
<li><p>fu = dot(N1, S) + d1</p></li>
<li><p>fv = dot(N2, S) + d2</p></li></ul>
<p>–</p>
<ul><li><p>err = fmax(fabs(fu), fabs(fv))</p></li>
<li><p>iterations += 1</p></li></ul>
<p>–</p>
<ul><li><p># t is the distance along the ray to this hit</p></li>
<li><p>cdef float t = distance(S, ray.org)</p></li></ul>
<p>+ cdef RayHitData hd = compute_patch_hit(patch.v, ray.org, ray.dir)</p>
<pre># only count this is it's the closest hit</pre>
<ul><li><p>if (t < ray.tnear or t > ray.Ng[0]):</p></li></ul>
<p>+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):</p>
<pre> return
</pre>
<ul><li><p>if (fabs(u) <= 1.0 and fabs(v) <= 1.0 and iterations < max_iter):</p></li></ul>
<p>+ if (fabs(hd.u) <= 1.0 and fabs(hd.v) <= 1.0 and hd.converged):</p>
<pre># we have a hit, so update ray information</pre>
<ul><li><p>ray.u = u</p></li>
<li><p>ray.v = v</p></li></ul>
<p>+ ray.u = hd.u + ray.v = hd.v</p>
<pre>ray.geomID = patch.geomID
ray.primID = item</pre>
<ul><li><p>ray.Ng[0] = t</p></li></ul>
<p>+ ray.Ng[0] = hd.t</p>
<pre># sample the solution at the calculated point
sample_hex20(patches, ray)</pre>
<p>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 @@</p>
<pre>from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
Patch</pre>
<p>-from yt.utilities.lib.mesh_intersection cimport patchSurfaceFunc +from yt.utilities.lib.primitives cimport patchSurfaceFunc</p>
<pre>from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \</pre>
<p>@@ -192,7 +192,7 @@</p>
<pre> elem_id = ray_id / 6
# fills "position" with the physical position of the hit</pre>
<ul><li><p>patchSurfaceFunc(data[ray_id], ray.u, ray.v, pos)</p></li></ul>
<p>+ patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)</p>
<pre> for i in range(3):
position[i] = <double> pos[i]
</pre>
<p>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 <strong>primitives, + const np.int64_t item, + np.float64_t[3] centroid) nogil + +cdef inline void triangle_bbox(const void <strong>primitives, + const np.int64_t item, + BBox</strong> 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</strong> primitives, + const np.int64_t item, + Ray* ray) nogil + +cdef inline void patch_centroid(const void <strong>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</strong> bbox) nogil</p>
<p>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) + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil: +# <a href="https://tavianator.com/fast-branchless-raybounding-box-intersections/">https://tavianator.com/fast-branchless-raybounding-box-intersections/</a> + + 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) + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline np.int64_t ray_triangle_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil: +# <a href="https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm">https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm</a> + + 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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void triangle_centroid(const void <strong>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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void triangle_bbox(const void <strong>primitives, + const np.int64_t item, + BBox</strong> 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]) + + +@cython.boundscheck(False) +@cython.wraparound(False) +@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)</strong>( 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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@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] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline np.int64_t ray_patch_intersect(const void* primitives, + const np.int64_t item, + Ray* ray) nogil: + + cdef 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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void patch_centroid(const void <strong>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 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline void patch_bbox(const void *primitives, + const np.int64_t item, + BBox</strong> 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])</p>
<p>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 @@</p>
<pre> from yt.utilities.lib import mesh_traversal
except ImportError:
mesh_traversal = NotAModule("pyembree")</pre>
<p>+ ytcfg["yt", “ray_tracing_engine”] = “yt”</p>
<pre>try:
from yt.utilities.lib import mesh_construction
except ImportError:
mesh_construction = NotAModule("pyembree")</pre>
<p>+ ytcfg["yt", “ray_tracing_engine”] = “yt”</p>
<pre>class RenderSource(ParallelAnalysisInterface):</pre>
<p>@@ -356,7 +358,7 @@</p>
<pre>self.field = field
self.volume = None
self.current_image = None</pre>
<ul><li><p>self.engine = ‘embree’</p></li></ul>
<p>+ self.engine = ytcfg.get("yt", “ray_tracing_engine”)</p>
<pre># default color map
self._cmap = ytcfg.get("yt", "default_colormap")</pre>
<p>@@ -374,11 +376,11 @@</p>
<pre> if self.engine == 'embree':
self.volume = mesh_traversal.YTEmbreeScene()
self.build_volume_embree()</pre>
<ul><li><p>elif self.engine == ‘bvh’:</p></li></ul>
<p>+ elif self.engine == ‘yt’:</p>
<pre>self.build_volume_bvh()
else:
raise NotImplementedError("Invalid ray-tracing engine selected. "</pre>
<ul><li><p>“Choices are ‘embree’ and 'bvh'.”)</p></li></ul>
<p>+ “Choices are ‘embree’ and 'yt'.”)</p>
<pre>def cmap():
'''</pre>
<p>@@ -491,15 +493,8 @@</p>
<pre> field_data = np.expand_dims(field_data, 1)
# Here, we decide whether to render based on high-order or</pre>
<ul><li><p># low-order geometry. Right now, high-order geometry is only</p></li>
<li><p># supported by Embree.</p></li>
<li><p>if indices.shape[1] == 20:</p></li>
<li><p># hexahedral</p></li>
<li><p>mylog.warning("20-node hexes not yet supported, " +</p></li>
<li><p>“dropping to 1st order.”)</p></li>
<li><p>field_data = field_data[:, 0:8]</p></li>
<li><p>indices = indices[:, 0:8]</p></li>
<li><p>elif indices.shape[1] == 27:</p></li></ul>
<p>+ # low-order geometry. + if indices.shape[1] == 27:</p>
<pre># hexahedral
mylog.warning("27-node hexes not yet supported, " +
"dropping to 1st order.")</pre>
<p>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 @@</p>
<pre>MeshSource, \
Scene, \
create_scene</pre>
<p>– – -@requires_module("pyembree") -def test_surface_mesh_render(): –</p>
<ul><li><p>images = []</p></li></ul>
<p>–</p>
<ul><li><p>ds = fake_tetrahedral_ds()</p></li>
<li><p>sc = Scene()</p></li>
<li><p>for field in ds.field_list:</p></li>
<li><p>sc.add_source(MeshSource(ds, field))</p></li>
<li><p>sc.add_camera()</p></li>
<li><p>im = sc.render()</p></li>
<li><p>images.append(im)</p></li></ul>
<p>–</p>
<ul><li><p>ds = fake_hexahedral_ds()</p></li>
<li><p>for field in ds.field_list:</p></li>
<li><p>sc.add_source(MeshSource(ds, field))</p></li>
<li><p>sc.add_camera()</p></li>
<li><p>im = sc.render()</p></li>
<li><p>images.append(im)</p></li></ul>
<p>–</p>
<ul><li><p>return images</p></li></ul>
<p>– +from yt.config import \ + ytcfg</p>
<pre>def compare(ds, im, test_prefix, decimals=12):
def mesh_render_image_func(filename_prefix):</pre>
<p>@@ -57,64 +35,132 @@</p>
<pre> test.prefix = test_prefix
return test
</pre>
<p>+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 + +@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() + +</p>
<pre>hex8 = "MOOSE_sample_data/out.e-s010"
hex8_fields = [('connect1', 'diffused'), ('connect2', 'convected')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(hex8)
@requires_module("pyembree")</pre>
<p>+def test_hex8_render_pyembree(): + for field in hex8_fields: + yield hex8_render("embree", field) + +@requires_ds(hex8)</p>
<pre>def test_hex8_render():
for field in hex8_fields:</pre>
<ul><li><p>ds = data_dir_load(hex8, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_hex8_%s_%s” % field)</p></li></ul>
<p>+ yield hex8_render("yt", field)</p>
<pre>tet4 = "MOOSE_sample_data/high_order_elems_tet4_refine_out.e"
tet4_fields = [("connect1", "u")]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(tet4)
@requires_module("pyembree")</pre>
<p>+def test_tet4_render_pyembree(): + for field in tet4_fields: + yield tet4_render("embree", field) + +@requires_ds(tet4)</p>
<pre>def test_tet4_render():
for field in tet4_fields:</pre>
<ul><li><p>ds = data_dir_load(tet4, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_tet4_%s_%s” % field)</p></li></ul>
<p>+ yield tet4_render("yt", field)</p>
<pre>hex20 = "MOOSE_sample_data/mps_out.e"
hex20_fields = [('connect2', 'temp')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(hex20)
@requires_module("pyembree")</pre>
<p>+def test_hex20_render_pyembree(): + for field in hex20_fields: + yield hex20_render("embree", field) + +@requires_ds(hex20)</p>
<pre>def test_hex20_render():
for field in hex20_fields:</pre>
<ul><li><p>ds = data_dir_load(hex20, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_hex20_%s_%s” % field)</p></li></ul>
<p>+ yield hex20_render("yt", field)</p>
<pre>wedge6 = "MOOSE_sample_data/wedge_out.e"
wedge6_fields = [('connect1', 'diffused')]
</pre>
<p>+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])) +</p>
<pre>@requires_ds(wedge6)
@requires_module("pyembree")</pre>
<p>+def test_wedge6_render_pyembree(): + for field in wedge6_fields: + yield wedge6_render("embree", field) + +@requires_ds(wedge6)</p>
<pre>def test_wedge6_render():
for field in wedge6_fields:</pre>
<ul><li><p>ds = data_dir_load(wedge6, kwargs={'step':-1})</p></li>
<li><p>sc = create_scene(ds, field)</p></li>
<li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “render_answers_wedge6_%s_%s” % field)</p></li></ul>
<p>+ yield wedge6_render("yt", field)</p>
<p>– -@requires_ds(hex8) -@requires_module("pyembree") -def test_perspective_mesh_render(): +def perspective_mesh_render(engine): + ytcfg["yt", “ray_tracing_engine”] = engine</p>
<pre>ds = data_dir_load(hex8)
sc = create_scene(ds, ("connect2", "diffused"))</pre>
<p>–</p>
<pre>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')</pre>
<p>@@ -122,12 +168,19 @@</p>
<pre>cam.set_position(cam_pos, north_vector)
cam.resolution = (800, 800)
im = sc.render()</pre>
<ul><li><p>yield compare(ds, im, “perspective_mesh_render”)</p></li></ul>
<p>– + return compare(ds, im, “%s_perspective_mesh_render” % engine) + +@requires_ds(hex8) +@requires_module("pyembree") +def test_perspective_mesh_render_pyembree(): + yield perspective_mesh_render("embree")</p>
<pre>@requires_ds(hex8)</pre>
<p>-@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</p>
<pre>ds = data_dir_load(hex8)
sc = Scene()
cam = sc.add_camera(ds)</pre>
<p>@@ -136,12 +189,18 @@</p>
<pre>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)</pre>
<p>–</p>
<pre>ms1 = MeshSource(ds, ('connect1', 'diffused'))
ms2 = MeshSource(ds, ('connect2', 'diffused'))</pre>
<p>–</p>
<pre>sc.add_source(ms1)
sc.add_source(ms2)</pre>
<p>+ im = sc.render() + return compare(ds, im, “%s_composite_mesh_render” % engine) + +@requires_ds(hex8) +@requires_module("pyembree") +def test_composite_mesh_render_pyembree(): + yield composite_mesh_render("embree")</p>
<ul><li><p>im = sc.render()</p></li>
<li><p>yield compare(ds, im, “composite_mesh_render”)</p></li></ul>
<p>+@requires_ds(hex8) +def test_composite_mesh_render(): + yield composite_mesh_render("yt")</p>
<p>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 @@</p>
<pre>kwargs = {'lens_type': params['lens_type']}
if engine == 'embree':
sampler = mesh_traversal.EmbreeMeshSampler(*args, **kwargs)</pre>
<ul><li><p>elif engine == ‘bvh’:</p></li></ul>
<p>+ elif engine == ‘yt’:</p>
<pre> sampler = mesh_traversal.BVHMeshSampler(*args, **kwargs)
return sampler</pre>
<p>Repository URL: <a href="https://bitbucket.org/yt_analysis/yt/">https://bitbucket.org/yt_analysis/yt/</a></p>
<p>—</p>
<p>This is a commit notification from bitbucket.org. You are receiving this because you have the service enabled, addressing the recipient of this email.</p>
<img src="http://link.bitbucket.org/wf/open?upn=ll4ctv0L-2ByeRZFC1LslHcg6aJmnQ70VruLbmeLQr27AEODQRKWrqz8ljagd-2F9PCW-2BDQbIqiweKYkWDf8FMN9cnl1c8pQDXgg2y2s2Qgrc3NJPtQwZi4RmCMcgNOtj-2BYujsz7sbdiwHPuma-2Fv8J-2BbPVHhkjWYaoyCL5dqBzHvopBwDkrojPhH14nrw4T5tT9jfYtr1klSM3rxdVqlmANQ9gw7D9BL4EKsgQSF0qygKY0-3D" alt="" width="1" height="1" border="0" style="height:1px !important;width:1px !important;border-width:0 !important;margin-top:0 !important;margin-bottom:0 !important;margin-right:0 !important;margin-left:0 !important;padding-top:0 !important;padding-bottom:0 !important;padding-right:0 !important;padding-left:0 !important;"/>
</body></html>