<html><body>
<p>1 new commit in yt:</p>
<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-2ByeRZFC1LslHcg6aJmnQ70VruLbmeLQr27CxTsX6zwy5ECxY5HFFrqYn-2FKZWA-2FYcrnNkre3v8D2GU3OzpN5o3Pr02d6kuZxnh529q4KwqZ1xJZWnC9hEsfS1-2FrzVRt0T2qJhrA-2BrV18liO7-2Bix4DjrwuH3Eo2zHEK39sOg58nxDZJY0Ez73jvT2KDxCrb8sHo7xzF-2B-2BuOhz9hCS7P1AkCtNNcvZLf1Izx3E-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>