[yt-svn] commit/yt: ngoldbaum: Merged in al007/yt (pull request #2401)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Oct 17 14:20:11 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/1d43d35f96cf/
Changeset: 1d43d35f96cf
Branch: yt
User: ngoldbaum
Date: 2016-10-17 21:19:44+00:00
Summary: Merged in al007/yt (pull request #2401)
Implement ability to volume render 3D meshes composed of second-order tetrahedrons
Affected #: 20 files
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -50,7 +50,7 @@
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_005:
+ local_varia_006:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
import numpy as np
cimport numpy as np
from yt.utilities.lib.element_mappings cimport ElementSampler
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
# node for the bounding volume hierarchy
cdef struct BVHNode:
@@ -69,8 +70,11 @@
cdef void _set_up_patches(self,
np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil
+ cdef void _set_up_tet_patches(self,
+ np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil
cdef void intersect(self, Ray* ray) nogil
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil
cdef void _recursive_intersect(self, Ray* ray, BVHNode* node) nogil
cdef BVHNode* _recursive_build(self, np.int64_t begin, np.int64_t end) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -18,19 +18,27 @@
Patch, \
ray_patch_intersect, \
patch_centroid, \
- patch_bbox
+ patch_bbox, \
+ TetPatch, \
+ ray_tet_patch_intersect, \
+ tet_patch_centroid, \
+ tet_patch_bbox
+
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
Q1Sampler3D, \
P1Sampler3D, \
W1Sampler3D, \
- S2Sampler3D
+ S2Sampler3D, \
+ Tet2Sampler3D
+
from yt.utilities.lib.vec3_ops cimport L2_norm
cdef ElementSampler Q1Sampler = Q1Sampler3D()
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
cdef extern from "platform_dep.h" nogil:
double fmax(double x, double y)
@@ -45,11 +53,11 @@
'''
This class implements a bounding volume hierarchy (BVH), a spatial acceleration
- structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
- partitioning the *volume* of the parent to create the children, we partition the
+ structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
+ partitioning the *volume* of the parent to create the children, we partition the
triangles themselves into 'left' or 'right' sub-trees. The bounding volume for a
node is then determined by computing the bounding volume of the triangles that
- belong to it. This allows us to quickly discard triangles that are not close
+ belong to it. This allows us to quickly discard triangles that are not close
to intersecting a given ray.
'''
@@ -82,6 +90,9 @@
elif self.num_verts_per_elem == 20:
self.num_prim_per_elem = 6
self.sampler = S2Sampler
+ elif self.num_verts_per_elem == 10:
+ self.num_prim_per_elem = 4
+ self.sampler = Tet2Sampler
else:
raise NotImplementedError("Could not determine element type for "
"nverts = %d. " % self.num_verts_per_elem)
@@ -109,7 +120,7 @@
self.vertices[vertex_offset + k] = vertices[indices[i,j]][k]
field_offset = i*self.num_field_per_elem
for j in range(self.num_field_per_elem):
- self.field_data[field_offset + j] = field_data[i][j]
+ self.field_data[field_offset + j] = field_data[i][j]
# set up primitives
if self.num_verts_per_elem == 20:
@@ -118,13 +129,19 @@
self.get_bbox = patch_bbox
self.get_intersect = ray_patch_intersect
self._set_up_patches(vertices, indices)
+ elif self.num_verts_per_elem == 10:
+ self.primitives = malloc(self.num_prim * sizeof(TetPatch))
+ self.get_centroid = tet_patch_centroid
+ self.get_bbox = tet_patch_bbox
+ self.get_intersect = ray_tet_patch_intersect
+ self._set_up_tet_patches(vertices, indices)
else:
self.primitives = malloc(self.num_prim * sizeof(Triangle))
self.get_centroid = triangle_centroid
self.get_bbox = triangle_bbox
self.get_intersect = ray_triangle_intersect
self._set_up_triangles(vertices, indices)
-
+
self.root = self._recursive_build(0, self.num_prim)
@cython.boundscheck(False)
@@ -156,6 +173,32 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+ cdef void _set_up_tet_patches(self, np.float64_t[:, :] vertices,
+ np.int64_t[:, :] indices) nogil:
+ cdef TetPatch* tet_patch
+ cdef np.int64_t i, j, k, ind, idim
+ cdef np.int64_t offset, prim_index
+ for i in range(self.num_elem):
+ offset = self.num_prim_per_elem*i
+ for j in range(self.num_prim_per_elem): # for each face
+ prim_index = offset + j
+ tet_patch = &( <TetPatch*> self.primitives)[prim_index]
+ self.prim_ids[prim_index] = prim_index
+ tet_patch.elem_id = i
+ for k in range(6): # for each vertex
+ ind = tet10_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ tet_patch.v[k][idim] = vertices[indices[i, ind]][idim]
+ self.get_centroid(self.primitives,
+ prim_index,
+ self.centroids[prim_index])
+ self.get_bbox(self.primitives,
+ prim_index,
+ &(self.bboxes[prim_index]))
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
cdef void _set_up_triangles(self, np.float64_t[:, :] vertices,
np.int64_t[:, :] indices) nogil:
# fill our array of primitives
@@ -181,7 +224,7 @@
tri_index,
self.centroids[tri_index])
self.get_bbox(self.primitives,
- tri_index,
+ tri_index,
&(self.bboxes[tri_index]))
cdef void _recursive_free(self, BVHNode* node) nogil:
@@ -191,6 +234,8 @@
free(node)
def __dealloc__(self):
+ if self.root == NULL:
+ return
self._recursive_free(self.root)
free(self.primitives)
free(self.prim_ids)
@@ -224,11 +269,11 @@
mid += 1
begin += 1
return mid
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void _get_node_bbox(self, BVHNode* node,
+ cdef void _get_node_bbox(self, BVHNode* node,
np.int64_t begin, np.int64_t end) nogil:
cdef np.int64_t i, j
cdef BBox box = self.bboxes[begin]
@@ -236,7 +281,7 @@
for j in range(3):
box.left_edge[j] = fmin(box.left_edge[j],
self.bboxes[i].left_edge[j])
- box.right_edge[j] = fmax(box.right_edge[j],
+ box.right_edge[j] = fmax(box.right_edge[j],
self.bboxes[i].right_edge[j])
node.bbox = box
@@ -245,7 +290,7 @@
@cython.cdivision(True)
cdef void intersect(self, Ray* ray) nogil:
self._recursive_intersect(ray, self.root)
-
+
if ray.elem_id < 0:
return
@@ -253,7 +298,7 @@
cdef np.int64_t i
for i in range(3):
position[i] = ray.origin[i] + ray.t_far*ray.direction[i]
-
+
cdef np.float64_t* vertex_ptr
cdef np.float64_t* field_ptr
vertex_ptr = self.vertices + ray.elem_id*self.num_verts_per_elem*3
@@ -297,17 +342,17 @@
node.end = end
self._get_node_bbox(node, begin, end)
-
+
# check for leaf
if (end - begin) <= LEAF_SIZE:
return node
-
+
# we use the "split in the middle of the longest axis approach"
# see: http://www.vadimkravcenko.com/bvh-tree-building/
# compute longest dimension
cdef np.int64_t ax = 0
- cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
+ cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
node.bbox.left_edge[0])
if fabs(node.bbox.right_edge[1] - node.bbox.left_edge[1]) > d:
ax = 1
@@ -323,7 +368,7 @@
if(mid == begin or mid == end):
mid = begin + (end-begin)/2
-
+
# recursively build sub-trees
node.left = self._recursive_build(begin, mid)
node.right = self._recursive_build(mid, end)
@@ -337,16 +382,16 @@
cdef void cast_rays(np.float64_t* image,
const np.float64_t* origins,
const np.float64_t* direction,
- const int N,
+ const int N,
BVH bvh) nogil:
- cdef Ray* ray
+ cdef Ray* ray
cdef int i, j, k
-
+
with nogil, parallel():
-
+
ray = <Ray *> malloc(sizeof(Ray))
-
+
for k in range(3):
ray.direction[k] = direction[k]
ray.inv_dir[k] = 1.0 / direction[k]
@@ -366,11 +411,11 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
+def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
np.ndarray[np.float64_t, ndim=2] origins,
np.ndarray[np.float64_t, ndim=1] direction,
BVH bvh):
-
+
cdef int N = origins.shape[0]
cast_rays(&image[0], &origins[0, 0], &direction[0], N, bvh)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -834,11 +834,10 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tris, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(2):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] - 1.0 > self.inclusion_tol):
+ return 0
return 1
cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
@@ -893,13 +892,39 @@
cdef int check_inside(self, double* mapped_coord) nogil:
# for canonical tets, we check whether the mapped_coords are between
# 0 and 1.
- cdef int i
- for i in range(3):
- if (mapped_coord[i] < -self.inclusion_tol or
- mapped_coord[i] - 1.0 > self.inclusion_tol):
- return 0
+ if (mapped_coord[0] < -self.inclusion_tol or \
+ mapped_coord[1] < -self.inclusion_tol or \
+ mapped_coord[2] < -self.inclusion_tol or \
+ mapped_coord[0] + mapped_coord[1] + mapped_coord[2] - 1.0 > \
+ self.inclusion_tol):
+ return 0
return 1
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_mesh_lines(self, double* mapped_coord) nogil:
+ cdef double u, v
+ cdef double thresh = 2.0e-2
+ if mapped_coord[0] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ elif mapped_coord[1] == 0:
+ u = mapped_coord[2]
+ v = mapped_coord[0]
+ elif mapped_coord[2] == 0:
+ u = mapped_coord[1]
+ v = mapped_coord[0]
+ else:
+ u = mapped_coord[1]
+ v = mapped_coord[2]
+ if ((u < thresh) or
+ (v < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh)):
+ return 1
+ return -1
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_construction.pxd
--- a/yt/utilities/lib/mesh_construction.pxd
+++ b/yt/utilities/lib/mesh_construction.pxd
@@ -18,3 +18,9 @@
unsigned int geomID
np.float64_t* vertices
np.float64_t* field_data
+
+ctypedef struct Tet_Patch:
+ float[6][3] v
+ unsigned int geomID
+ np.float64_t* vertices
+ np.float64_t* field_data
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -1,5 +1,5 @@
"""
-This file contains the ElementMesh classes, which represent the target that the
+This file contains the ElementMesh classes, which represent the target that the
rays will be cast at when rendering finite element data. This class handles
the interface between the internal representation of the mesh and the pyembree
representation.
@@ -21,7 +21,7 @@
from libc.math cimport fmax, sqrt
cimport numpy as np
-cimport pyembree.rtcore as rtc
+cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry_user as rtcgu
@@ -37,13 +37,15 @@
sample_wedge
from mesh_intersection cimport \
patchIntersectFunc, \
- patchBoundsFunc
+ patchBoundsFunc, \
+ tet_patchIntersectFunc, \
+ tet_patchBoundsFunc
from yt.utilities.exceptions import YTElementTypeNotRecognized
cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
-
+
int HEX_NV
int HEX_NT
int TETRA_NV
@@ -54,13 +56,14 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef class LinearElementMesh:
r'''
This creates a 1st-order mesh to be ray-traced with embree.
- Currently, we handle non-triangular mesh types by converting them
+ Currently, we handle non-triangular mesh types by converting them
to triangular meshes. This class performs this transformation.
Currently, this is implemented for hexahedral and tetrahedral
meshes.
@@ -71,21 +74,21 @@
scene : EmbreeScene
This is the scene to which the constructed polygons will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the polygon mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the polygon mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
- This should either have the shape (num_elements, 4) or
- (num_elements, 8) for tetrahedral and hexahedral meshes,
- respectively. For tetrahedral meshes, each element will
+ This should either have the shape (num_elements, 4) or
+ (num_elements, 8) for tetrahedral and hexahedral meshes,
+ respectively. For tetrahedral meshes, each element will
be represented by four triangles in the scene. For hex meshes,
- each element will be represented by 12 triangles, 2 for each
+ each element will be represented by 12 triangles, 2 for each
face. For hex meshes, we assume that the node ordering is as
- defined here:
+ defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
cdef Vertex* vertices
@@ -93,7 +96,7 @@
cdef unsigned int mesh
cdef double* field_data
cdef rtcg.RTCFilterFunc filter_func
- # triangles per element, vertices per element, and field points per
+ # triangles per element, vertices per element, and field points per
# element, respectively:
cdef int tpe, vpe, fpe
cdef int[MAX_NUM_TRI][3] tri_array
@@ -101,7 +104,7 @@
cdef MeshDataContainer datac
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray data):
@@ -119,7 +122,7 @@
self.tpe = TETRA_NT
self.tri_array = triangulate_tetra
else:
- raise YTElementTypeNotRecognized(vertices.shape[1],
+ raise YTElementTypeNotRecognized(vertices.shape[1],
indices.shape[1])
self._build_from_indices(scene, vertices, indices)
@@ -142,7 +145,7 @@
for i in range(nv):
vertices[i].x = vertices_in[i, 0]
vertices[i].y = vertices_in[i, 1]
- vertices[i].z = vertices_in[i, 2]
+ vertices[i].z = vertices_in[i, 2]
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_VERTEX_BUFFER,
vertices, 0, sizeof(Vertex))
@@ -156,7 +159,7 @@
rtcg.rtcSetBuffer(scene.scene_i, mesh, rtcg.RTC_INDEX_BUFFER,
triangles, 0, sizeof(Triangle))
- cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
+ cdef int* element_indices = <int *> malloc(ne * self.vpe * sizeof(int))
for i in range(ne):
for j in range(self.vpe):
element_indices[i*self.vpe + j] = indices_in[i][j]
@@ -189,7 +192,7 @@
datac.vpe = self.vpe
datac.fpe = self.fpe
self.datac = datac
-
+
rtcg.rtcSetUserData(scene.scene_i, self.mesh, &self.datac)
cdef void _set_sampler_type(self, YTEmbreeScene scene):
@@ -206,7 +209,7 @@
rtcg.rtcSetIntersectionFilterFunction(scene.scene_i,
self.mesh,
self.filter_func)
-
+
def __dealloc__(self):
free(self.field_data)
free(self.element_indices)
@@ -227,84 +230,137 @@
scene : EmbreeScene
This is the scene to which the constructed patches will be
added.
- vertices : a np.ndarray of floats.
- This specifies the x, y, and z coordinates of the vertices in
- the mesh. This should either have the shape
- (num_vertices, 3). For example, vertices[2][1] should give the
+ vertices : a np.ndarray of floats.
+ This specifies the x, y, and z coordinates of the vertices in
+ the mesh. This should either have the shape
+ (num_vertices, 3). For example, vertices[2][1] should give the
y-coordinate of the 3rd vertex in the mesh.
indices : a np.ndarray of ints
This should have the shape (num_elements, 20). Each hex will be
- represented in the scene by 6 bi-quadratic patches. We assume that
- the node ordering is as defined here:
+ represented in the scene by 6 bi-quadratic patches. We assume that
+ the node ordering is as defined here:
http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
-
+
'''
- cdef Patch* patches
+ cdef void* patches
cdef np.float64_t* vertices
cdef np.float64_t* field_data
cdef unsigned int mesh
- # patches per element, vertices per element, and field points per
- # element, respectively:
- cdef int ppe, vpe, fpe
+ # patches per element, vertices per element, vertices per face,
+ # and field points per element, respectively:
+ cdef int ppe, vpe, vpf, fpe
def __init__(self, YTEmbreeScene scene,
- np.ndarray vertices,
+ np.ndarray vertices,
np.ndarray indices,
np.ndarray field_data):
+ cdef int i, j
- # only 20-point hexes are supported right now.
+ # 20-point hexes
if indices.shape[1] == 20:
self.vpe = 20
+ self.ppe = 6
+ self.vpf = 8
+ self._build_from_indices_hex20(scene, vertices, indices, field_data)
+ # 10-point tets
+ elif indices.shape[1] == 10:
+ self.vpe = 10
+ self.ppe = 4
+ self.vpf = 6
+ self._build_from_indices_tet10(scene, vertices, indices, field_data)
else:
raise NotImplementedError
- self._build_from_indices(scene, vertices, indices, field_data)
-
- cdef void _build_from_indices(self, YTEmbreeScene scene,
+ cdef void _build_from_indices_hex20(self, YTEmbreeScene scene,
np.ndarray vertices_in,
np.ndarray indices_in,
np.ndarray field_data):
- cdef int i, j, ind, idim
+ cdef int i, j, k, ind, idim
cdef int ne = indices_in.shape[0]
cdef int nv = vertices_in.shape[0]
- cdef int npatch = 6*ne;
+ cdef int npatch = self.ppe*ne;
cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
cdef np.ndarray[np.float64_t, ndim=2] element_vertices
cdef Patch* patches = <Patch*> malloc(npatch * sizeof(Patch))
- self.vertices = <np.float64_t*> malloc(20 * ne * 3 * sizeof(np.float64_t))
- self.field_data = <np.float64_t*> malloc(20 * ne * sizeof(np.float64_t))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
for i in range(ne):
element_vertices = vertices_in[indices_in[i]]
- for j in range(20):
- self.field_data[i*20 + j] = field_data[i][j]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
for k in range(3):
- self.vertices[i*20*3 + j*3 + k] = element_vertices[j][k]
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
cdef Patch* patch
for i in range(ne): # for each element
element_vertices = vertices_in[indices_in[i]]
- for j in range(6): # for each face
- patch = &(patches[i*6+j])
+ for j in range(self.ppe): # for each face
+ patch = &(patches[i*self.ppe+j])
patch.geomID = mesh
- for k in range(8): # for each vertex
+ for k in range(self.vpf): # for each vertex
ind = hex20_faces[j][k]
for idim in range(3): # for each spatial dimension (yikes)
patch.v[k][idim] = element_vertices[ind][idim]
- patch.vertices = self.vertices + i*20*3
- patch.field_data = self.field_data + i*20
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
self.patches = patches
self.mesh = mesh
- rtcg.rtcSetUserData(scene.scene_i, self.mesh, self.patches)
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
<rtcgu.RTCBoundsFunc> patchBoundsFunc)
rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
<rtcgu.RTCIntersectFunc> patchIntersectFunc)
+ cdef void _build_from_indices_tet10(self, YTEmbreeScene scene,
+ np.ndarray vertices_in,
+ np.ndarray indices_in,
+ np.ndarray field_data):
+ cdef int i, j, k, ind, idim
+ cdef int ne = indices_in.shape[0]
+ cdef int nv = vertices_in.shape[0]
+ cdef int npatch = self.ppe*ne;
+
+ cdef unsigned int mesh = rtcgu.rtcNewUserGeometry(scene.scene_i, npatch)
+ cdef np.ndarray[np.float64_t, ndim=2] element_vertices
+ cdef Tet_Patch* patches = <Tet_Patch*> malloc(npatch * sizeof(Tet_Patch))
+ self.vertices = <np.float64_t*> malloc(self.vpe * ne * 3 * sizeof(np.float64_t))
+ self.field_data = <np.float64_t*> malloc(self.vpe * ne * sizeof(np.float64_t))
+
+ for i in range(ne):
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.vpe):
+ self.field_data[i*self.vpe + j] = field_data[i][j]
+ for k in range(3):
+ self.vertices[i*self.vpe*3 + j*3 + k] = element_vertices[j][k]
+
+ cdef Tet_Patch* patch
+ for i in range(ne): # for each element
+ element_vertices = vertices_in[indices_in[i]]
+ for j in range(self.ppe): # for each face
+ patch = &(patches[i*self.ppe+j])
+ patch.geomID = mesh
+ for k in range(self.vpf): # for each vertex
+ ind = tet10_faces[j][k]
+ for idim in range(3): # for each spatial dimension (yikes)
+ patch.v[k][idim] = element_vertices[ind][idim]
+ patch.vertices = self.vertices + i*self.vpe*3
+ patch.field_data = self.field_data + i*self.vpe
+
+ self.patches = patches
+ self.mesh = mesh
+
+ rtcg.rtcSetUserData(scene.scene_i, self.mesh, patches)
+ rtcgu.rtcSetBoundsFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCBoundsFunc> tet_patchBoundsFunc)
+ rtcgu.rtcSetIntersectFunction(scene.scene_i, self.mesh,
+ <rtcgu.RTCIntersectFunc> tet_patchIntersectFunc)
+
+
def __dealloc__(self):
free(self.patches)
free(self.vertices)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_intersection.pxd
--- a/yt/utilities/lib/mesh_intersection.pxd
+++ b/yt/utilities/lib/mesh_intersection.pxd
@@ -1,7 +1,7 @@
cimport pyembree.rtcore as rtc
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry as rtcg
-from yt.utilities.lib.mesh_construction cimport Patch
+from yt.utilities.lib.mesh_construction cimport Patch, Tet_Patch
cimport cython
cdef void patchIntersectFunc(Patch* patches,
@@ -11,3 +11,11 @@
cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil
+
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil
+
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_intersection.pyx
--- a/yt/utilities/lib/mesh_intersection.pyx
+++ b/yt/utilities/lib/mesh_intersection.pyx
@@ -20,26 +20,30 @@
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmin, fmax, sqrt
-from yt.utilities.lib.mesh_samplers cimport sample_hex20
+from yt.utilities.lib.mesh_samplers cimport sample_hex20, sample_tet10
from yt.utilities.lib.bounding_volume_hierarchy cimport BBox
from yt.utilities.lib.primitives cimport \
patchSurfaceFunc, \
patchSurfaceDerivU, \
patchSurfaceDerivV, \
RayHitData, \
- compute_patch_hit
+ compute_patch_hit, \
+ tet_patchSurfaceFunc, \
+ tet_patchSurfaceDerivU, \
+ tet_patchSurfaceDerivV, \
+ compute_tet_patch_hit
from vec3_ops cimport dot, subtract, cross, distance
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void patchBoundsFunc(Patch* patches,
+cdef void patchBoundsFunc(Patch* patches,
size_t item,
rtcg.RTCBounds* bounds_o) nogil:
cdef Patch patch = patches[item]
-
+
cdef float lo_x = 1.0e300
cdef float lo_y = 1.0e300
cdef float lo_z = 1.0e300
@@ -88,8 +92,71 @@
ray.geomID = patch.geomID
ray.primID = item
ray.Ng[0] = hd.t
-
+
# sample the solution at the calculated point
sample_hex20(patches, ray)
return
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchBoundsFunc(Tet_Patch* tet_patches,
+ size_t item,
+ rtcg.RTCBounds* bounds_o) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef float lo_x = 1.0e300
+ cdef float lo_y = 1.0e300
+ cdef float lo_z = 1.0e300
+
+ cdef float hi_x = -1.0e300
+ cdef float hi_y = -1.0e300
+ cdef float hi_z = -1.0e300
+
+ cdef int i
+ for i in range(6):
+ lo_x = fmin(lo_x, tet_patch.v[i][0])
+ lo_y = fmin(lo_y, tet_patch.v[i][1])
+ lo_z = fmin(lo_z, tet_patch.v[i][2])
+ hi_x = fmax(hi_x, tet_patch.v[i][0])
+ hi_y = fmax(hi_y, tet_patch.v[i][1])
+ hi_z = fmax(hi_z, tet_patch.v[i][2])
+
+ bounds_o.lower_x = lo_x
+ bounds_o.lower_y = lo_y
+ bounds_o.lower_z = lo_z
+ bounds_o.upper_x = hi_x
+ bounds_o.upper_y = hi_y
+ bounds_o.upper_z = hi_z
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchIntersectFunc(Tet_Patch* tet_patches,
+ rtcr.RTCRay& ray,
+ size_t item) nogil:
+
+ cdef Tet_Patch tet_patch = tet_patches[item]
+
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.org, ray.dir)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.tnear or hd.t > ray.Ng[0]):
+ return
+
+ if (hd.converged and 0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1):
+
+ # we have a hit, so update ray information
+ ray.u = hd.u
+ ray.v = hd.v
+ ray.geomID = tet_patch.geomID
+ ray.primID = item
+ ray.Ng[0] = hd.t
+
+ # sample the solution at the calculated point
+ sample_tet10(tet_patches, ray)
+
+ return
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_samplers.pxd
--- a/yt/utilities/lib/mesh_samplers.pxd
+++ b/yt/utilities/lib/mesh_samplers.pxd
@@ -13,3 +13,6 @@
cdef void sample_hex20(void* userPtr,
rtcr.RTCRay& ray) nogil
+
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_samplers.pyx
--- a/yt/utilities/lib/mesh_samplers.pyx
+++ b/yt/utilities/lib/mesh_samplers.pyx
@@ -18,14 +18,16 @@
from pyembree.rtcore cimport Vec3f, Triangle, Vertex
from yt.utilities.lib.mesh_construction cimport \
MeshDataContainer, \
- Patch
-from yt.utilities.lib.primitives cimport patchSurfaceFunc
+ Patch, \
+ Tet_Patch
+from yt.utilities.lib.primitives cimport patchSurfaceFunc, tet_patchSurfaceFunc
from yt.utilities.lib.element_mappings cimport \
ElementSampler, \
P1Sampler3D, \
Q1Sampler3D, \
S2Sampler3D, \
- W1Sampler3D
+ W1Sampler3D, \
+ Tet2Sampler3D
cimport numpy as np
cimport cython
from libc.math cimport fabs, fmax
@@ -34,6 +36,7 @@
cdef ElementSampler P1Sampler = P1Sampler3D()
cdef ElementSampler S2Sampler = S2Sampler3D()
cdef ElementSampler W1Sampler = W1Sampler3D()
+cdef ElementSampler Tet2Sampler = Tet2Sampler3D()
@cython.boundscheck(False)
@@ -94,7 +97,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(8):
element_indices[i] = data.element_indices[elem_id*8+i]
@@ -104,7 +107,7 @@
for i in range(8):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -145,7 +148,7 @@
elem_id = ray_id / data.tpe
get_hit_position(position, userPtr, ray)
-
+
for i in range(6):
element_indices[i] = data.element_indices[elem_id*6+i]
@@ -155,7 +158,7 @@
for i in range(6):
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
@@ -195,14 +198,14 @@
patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
for i in range(3):
position[i] = <double> pos[i]
-
+
# we use ray.time to pass the value of the field
cdef double mapped_coord[3]
S2Sampler.map_real_to_unit(mapped_coord, patch.vertices, position)
val = S2Sampler.sample_at_unit_point(mapped_coord, patch.field_data)
ray.time = val
ray.instID = S2Sampler.check_mesh_lines(mapped_coord)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -228,14 +231,14 @@
# ray_id records the id number of the hit according to
# embree, in which the primitives are triangles. Here,
# we convert this to the element id by dividing by the
- # number of triangles per element.
+ # number of triangles per element.
elem_id = ray_id / data.tpe
for i in range(4):
element_indices[i] = data.element_indices[elem_id*4+i]
vertices[i*3] = data.vertices[element_indices[i]].x
vertices[i*3 + 1] = data.vertices[element_indices[i]].y
- vertices[i*3 + 2] = data.vertices[element_indices[i]].z
+ vertices[i*3 + 2] = data.vertices[element_indices[i]].z
for i in range(data.fpe):
field_data[i] = data.field_data[elem_id*data.fpe+i]
@@ -249,3 +252,39 @@
val = P1Sampler.sample_at_unit_point(mapped_coord, field_data)
ray.time = val
ray.instID = P1Sampler.check_mesh_lines(mapped_coord)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.initializedcheck(False)
+ at cython.cdivision(True)
+cdef void sample_tet10(void* userPtr,
+ rtcr.RTCRay& ray) nogil:
+ cdef int ray_id, elem_id, i
+ cdef double val
+ cdef double[3] position
+ cdef float[3] pos
+ cdef Tet_Patch* data
+
+ data = <Tet_Patch*> userPtr
+ ray_id = ray.primID
+ if ray_id == -1:
+ return
+ cdef Tet_Patch tet_patch = data[ray_id]
+
+ # ray_id records the id number of the hit according to
+ # embree, in which the primitives are patches. Here,
+ # we convert this to the element id by dividing by the
+ # number of patches per element.
+ elem_id = ray_id / 4
+
+ # fills "position" with the physical position of the hit
+ tet_patchSurfaceFunc(data[ray_id].v, ray.u, ray.v, pos)
+ for i in range(3):
+ position[i] = <double> pos[i]
+
+ # we use ray.time to pass the value of the field
+ cdef double mapped_coord[3]
+ Tet2Sampler.map_real_to_unit(mapped_coord, tet_patch.vertices, position)
+ val = Tet2Sampler.sample_at_unit_point(mapped_coord, tet_patch.field_data)
+ ray.time = val
+ ray.instID = Tet2Sampler.check_mesh_lines(mapped_coord)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_triangulation.h
--- a/yt/utilities/lib/mesh_triangulation.h
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -12,7 +12,7 @@
// here: http://homepages.cae.wisc.edu/~tautges/papers/cnmev3.pdf
// Note that this is the case for Exodus II data.
int triangulate_hex[MAX_NUM_TRI][3] = {
- {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
+ {0, 2, 1}, {0, 3, 2}, // Face is 3 2 1 0
{4, 5, 6}, {4, 6, 7}, // Face is 4 5 6 7
{0, 1, 5}, {0, 5, 4}, // Face is 0 1 5 4
{1, 2, 6}, {1, 6, 5}, // Face is 1 2 6 5
@@ -22,7 +22,7 @@
// Similarly, this is used to triangulate the tetrahedral cells
int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 3},
+ {0, 1, 3},
{2, 3, 1},
{0, 3, 2},
{0, 2, 1},
@@ -57,10 +57,18 @@
// This is used to select faces from a 20-sided hex element
int hex20_faces[6][8] = {
- {0, 1, 5, 4, 12, 8, 13, 16},
+ {0, 1, 5, 4, 12, 8, 13, 16},
{1, 2, 6, 5, 13, 9, 14, 17},
{3, 2, 6, 7, 15, 10, 14, 18},
{0, 3, 7, 4, 12, 11, 15, 19},
{4, 5, 6, 7, 19, 16, 17, 18},
{0, 1, 2, 3, 11, 8, 9, 10}
};
+
+// This is used to select faces from a second-order tet element
+int tet10_faces[4][6] = {
+ {0, 1, 3, 4, 8, 7},
+ {2, 3, 1, 9, 8, 5},
+ {0, 3, 2, 7, 9, 6},
+ {0, 2, 1, 6, 5, 4}
+};
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/mesh_triangulation.pyx
--- a/yt/utilities/lib/mesh_triangulation.pyx
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -18,6 +18,7 @@
int triangulate_tetra[MAX_NUM_TRI][3]
int triangulate_wedge[MAX_NUM_TRI][3]
int hex20_faces[6][8]
+ int tet10_faces[4][6]
cdef struct TriNode:
np.uint64_t key
@@ -35,7 +36,7 @@
if not found:
return 0
return 1
-
+
cdef np.uint64_t triangle_hash(np.int64_t tri[3]) nogil:
# http://stackoverflow.com/questions/1536393/good-hash-function-for-permutations
cdef np.uint64_t h = 1
@@ -58,13 +59,13 @@
cdef TriNode **table
cdef np.uint64_t num_items
-
+
def __cinit__(self):
self.table = <TriNode**> malloc(TABLE_SIZE * sizeof(TriNode*))
for i in range(TABLE_SIZE):
self.table[i] = NULL
self.num_items = 0
-
+
def __dealloc__(self):
cdef np.int64_t i
cdef TriNode *node
@@ -77,7 +78,7 @@
free(delete_node)
self.table[i] = NULL
free(self.table)
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_exterior_tris(self):
@@ -102,7 +103,7 @@
element_map[counter] = node.elem
counter += 1
node = node.next_node
-
+
return tri_indices, element_map
cdef TriNode* _allocate_new_node(self,
@@ -118,13 +119,13 @@
new_node.next_node = NULL
self.num_items += 1
return new_node
-
+
@cython.cdivision(True)
cdef void update(self, np.int64_t tri[3], np.int64_t elem) nogil:
cdef np.uint64_t key = triangle_hash(tri)
cdef np.uint64_t index = key % TABLE_SIZE
cdef TriNode *node = self.table[index]
-
+
if node == NULL:
self.table[index] = self._allocate_new_node(tri, key, elem)
return
@@ -139,7 +140,7 @@
elif node.next_node == NULL:
node.next_node = self._allocate_new_node(tri, key, elem)
return
-
+
# walk through node list
cdef TriNode* prev = node
node = node.next_node
@@ -165,7 +166,7 @@
cdef np.int64_t VPE # num verts per element
cdef np.int64_t TPE # num tris per element
cdef int[MAX_NUM_TRI][3] tri_array
-
+
def __cinit__(self, np.int64_t[:, ::1] indices):
'''
@@ -179,7 +180,7 @@
if (self.VPE == 8 or self.VPE == 20 or self.VPE == 27):
self.TPE = HEX_NT
self.tri_array = triangulate_hex
- elif self.VPE == 4:
+ elif (self.VPE == 4 or self.VPE == 10):
self.TPE = TETRA_NT
self.tri_array = triangulate_tetra
elif self.VPE == 6:
@@ -190,7 +191,7 @@
self.num_tri = self.TPE * self.num_elem
self.num_verts = self.num_tri * 3
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def cull_interior_triangles(np.int64_t[:, ::1] indices):
@@ -213,7 +214,7 @@
s.update(tri, i)
return s.get_exterior_tris()
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
def get_vertex_data(np.float64_t[:, ::1] coords,
@@ -230,7 +231,7 @@
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t num_verts = coords.shape[0]
cdef np.float32_t[:] vertex_data = np.zeros(num_verts, dtype="float32")
-
+
cdef np.int64_t i, j
for i in range(m.num_elem):
for j in range(m.VPE):
@@ -256,17 +257,17 @@
cdef np.int64_t num_tri = exterior_tris.shape[0]
cdef np.int64_t num_verts = 3 * num_tri
cdef np.int64_t num_coords = 3 * num_verts
-
+
cdef np.float32_t[:] vertex_data
if data.ndim == 2:
vertex_data = get_vertex_data(coords, data, indices)
else:
vertex_data = data.astype("float32")
-
+
cdef np.int32_t[:] tri_indices = np.empty(num_verts, dtype=np.int32)
cdef np.float32_t[:] tri_data = np.empty(num_verts, dtype=np.float32)
cdef np.float32_t[:] tri_coords = np.empty(num_coords, dtype=np.float32)
-
+
cdef np.int64_t vert_index, i, j, k
for i in range(num_tri):
for j in range(3):
@@ -278,7 +279,7 @@
tri_indices[vert_index] = vert_index
for k in range(3):
tri_coords[vert_index*3 + k] = coords[exterior_tris[i, j], k]
-
+
return np.array(tri_coords), np.array(tri_data), np.array(tri_indices)
@cython.boundscheck(False)
@@ -291,10 +292,10 @@
coordinates and the data values.
'''
-
+
cdef MeshInfoHolder m = MeshInfoHolder(indices)
cdef np.int64_t[:, ::1] tri_indices = np.empty((m.num_tri, 3), dtype=np.int64)
-
+
cdef np.int64_t i, j, k
for i in range(m.num_elem):
for j in range(m.TPE):
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/primitives.pxd
--- a/yt/utilities/lib/primitives.pxd
+++ b/yt/utilities/lib/primitives.pxd
@@ -1,4 +1,4 @@
-cimport cython
+cimport cython
cimport cython.floating
import numpy as np
cimport numpy as np
@@ -66,7 +66,7 @@
cdef RayHitData compute_patch_hit(cython.floating[8][3] verts,
cython.floating[3] ray_origin,
cython.floating[3] ray_direction) nogil
-
+
cdef np.int64_t ray_patch_intersect(const void* primitives,
const np.int64_t item,
Ray* ray) nogil
@@ -78,3 +78,38 @@
cdef void patch_bbox(const void *primitives,
const np.int64_t item,
BBox* bbox) nogil
+
+cdef struct TetPatch:
+ np.float64_t[6][3] v # 6 vertices per patch
+ np.int64_t elem_id
+
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil
+
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil
+
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil
+
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil
+
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil
+
+cdef void tet_patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil
+
+cdef void tet_patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/lib/primitives.pyx
--- a/yt/utilities/lib/primitives.pyx
+++ b/yt/utilities/lib/primitives.pyx
@@ -21,15 +21,15 @@
cdef np.float64_t tmin = -INF
cdef np.float64_t tmax = INF
-
+
cdef np.int64_t i
cdef np.float64_t t1, t2
for i in range(3):
t1 = (bbox.left_edge[i] - ray.origin[i])*ray.inv_dir[i]
- t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
+ t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
tmin = fmax(tmin, fmin(t1, t2))
tmax = fmin(tmax, fmax(t1, t2))
-
+
return tmax >= fmax(tmin, 0.0)
@cython.boundscheck(False)
@@ -53,7 +53,7 @@
cdef np.float64_t det, inv_det
det = dot(e1, P)
- if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
+ if(det > -DETERMINANT_EPS and det < DETERMINANT_EPS):
return False
inv_det = 1.0 / det
@@ -87,7 +87,7 @@
cdef void triangle_centroid(const void *primitives,
const np.int64_t item,
np.float64_t[3] centroid) nogil:
-
+
cdef Triangle tri = (<Triangle*> primitives)[item]
cdef np.int64_t i
for i in range(3):
@@ -112,7 +112,7 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceFunc(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] S) nogil:
@@ -132,9 +132,9 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivU(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
- cython.floating[3] Su) nogil:
+ cython.floating[3] Su) nogil:
cdef int i
for i in range(3):
Su[i] = (-0.25*(v - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
@@ -149,10 +149,10 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void patchSurfaceDerivV(const cython.floating[8][3] verts,
- const cython.floating u,
+ const cython.floating u,
const cython.floating v,
cython.floating[3] Sv) nogil:
- cdef int i
+ cdef int i
for i in range(3):
Sv[i] = (-0.25*(u - 1.0)*(u + v + 1) - 0.25*(u - 1.0)*(v - 1.0))*verts[0][i] + \
(-0.25*(u + 1.0)*(u - v - 1) + 0.25*(u + 1.0)*(v - 1.0))*verts[1][i] + \
@@ -196,7 +196,7 @@
cdef cython.floating fu = dot(N1, S) + d1
cdef cython.floating fv = dot(N2, S) + d2
cdef cython.floating err = fmax(fabs(fu), fabs(fv))
-
+
# begin Newton interation
cdef cython.floating tol = 1.0e-5
cdef int iterations = 0
@@ -213,11 +213,11 @@
J21 = dot(N2, Su)
J22 = dot(N2, Sv)
det = (J11*J22 - J12*J21)
-
+
# update the u, v values
u -= ( J22*fu - J12*fv) / det
v -= (-J21*fu + J11*fv) / det
-
+
patchSurfaceFunc(verts, u, v, S)
fu = dot(N1, S) + d1
fv = dot(N2, S) + d2
@@ -227,7 +227,7 @@
# t is the distance along the ray to this hit
cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
-
+
# return hit data
cdef RayHitData hd
hd.u = u
@@ -247,7 +247,7 @@
cdef Patch patch = (<Patch*> primitives)[item]
cdef RayHitData hd = compute_patch_hit(patch.v, ray.origin, ray.direction)
-
+
# only count this is it's the closest hit
if (hd.t < ray.t_near or hd.t > ray.t_far):
return False
@@ -259,7 +259,7 @@
return True
return False
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -291,7 +291,7 @@
cdef np.int64_t i, j
cdef Patch patch = (<Patch*> primitives)[item]
-
+
for j in range(3):
bbox.left_edge[j] = patch.v[0][j]
bbox.right_edge[j] = patch.v[0][j]
@@ -300,3 +300,196 @@
for j in range(3):
bbox.left_edge[j] = fmin(bbox.left_edge[j], patch.v[i][j])
bbox.right_edge[j] = fmax(bbox.right_edge[j], patch.v[i][j])
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceFunc(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] S) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ S[i] = (1.0 - 3.0*u + 2.0*u*u - 3.0*v + 2.0*v*v + 4.0*u*v)*verts[0][i] + \
+ (-u + 2.0*u*u)*verts[1][i] + \
+ (-v + 2.0*v*v)*verts[2][i] + \
+ (4.0*u - 4.0*u*u - 4.0*u*v)*verts[3][i] + \
+ (4.0*u*v)*verts[4][i] + \
+ (4.0*v - 4.0*v*v - 4.0*u*v)*verts[5][i]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceDerivU(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Su) nogil:
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Su[i] = (-3.0 + 4.0*u + 4.0*v)*verts[0][i] + \
+ (-1.0 + 4.0*u)*verts[1][i] + \
+ (4.0 - 8.0*u - 4.0*v)*verts[3][i] + \
+ (4.0*v)*verts[4][i] + \
+ (-4.0*v)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patchSurfaceDerivV(const cython.floating[6][3] verts,
+ const cython.floating u,
+ const cython.floating v,
+ cython.floating[3] Sv) nogil:
+
+ cdef int i
+ # Computes for canonical triangle coordinates
+ for i in range(3):
+ Sv[i] = (-3.0 + 4.0*v + 4.0*u)*verts[0][i] + \
+ (-1.0 + 4.0*v)*verts[2][i] + \
+ (-4.0*u)*verts[3][i] + \
+ (4.0*u)*verts[4][i] + \
+ (4.0 - 8.0*v - 4.0*u)*verts[5][i]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef RayHitData compute_tet_patch_hit(cython.floating[6][3] verts,
+ cython.floating[3] ray_origin,
+ cython.floating[3] ray_direction) nogil:
+
+ # first we compute the two planes that define the ray.
+ cdef cython.floating[3] n, N1, N2
+ cdef cython.floating A = dot(ray_direction, ray_direction)
+ for i in range(3):
+ n[i] = ray_direction[i] / A
+
+ if ((fabs(n[0]) > fabs(n[1])) and (fabs(n[0]) > fabs(n[2]))):
+ N1[0] = n[1]
+ N1[1] =-n[0]
+ N1[2] = 0.0
+ else:
+ N1[0] = 0.0
+ N1[1] = n[2]
+ N1[2] =-n[1]
+ cross(N1, n, N2)
+
+ cdef cython.floating d1 = -dot(N1, ray_origin)
+ cdef cython.floating d2 = -dot(N2, ray_origin)
+
+ # the initial guess is set to zero
+ cdef cython.floating u = 0.0
+ cdef cython.floating v = 0.0
+ cdef cython.floating[3] S
+ tet_patchSurfaceFunc(verts, u, v, S)
+ cdef cython.floating fu = dot(N1, S) + d1
+ cdef cython.floating fv = dot(N2, S) + d2
+ cdef cython.floating err = fmax(fabs(fu), fabs(fv))
+
+ # begin Newton interation
+ cdef cython.floating tol = 1.0e-5
+ cdef int iterations = 0
+ cdef int max_iter = 10
+ cdef cython.floating[3] Su
+ cdef cython.floating[3] Sv
+ cdef cython.floating J11, J12, J21, J22, det
+ while ((err > tol) and (iterations < max_iter)):
+ # compute the Jacobian
+ tet_patchSurfaceDerivU(verts, u, v, Su)
+ tet_patchSurfaceDerivV(verts, u, v, Sv)
+ J11 = dot(N1, Su)
+ J12 = dot(N1, Sv)
+ J21 = dot(N2, Su)
+ J22 = dot(N2, Sv)
+ det = (J11*J22 - J12*J21)
+
+ # update the u, v values
+ u -= ( J22*fu - J12*fv) / det
+ v -= (-J21*fu + J11*fv) / det
+
+ tet_patchSurfaceFunc(verts, u, v, S)
+ fu = dot(N1, S) + d1
+ fv = dot(N2, S) + d2
+
+ err = fmax(fabs(fu), fabs(fv))
+ iterations += 1
+
+ # t is the distance along the ray to this hit
+ cdef cython.floating t = distance(S, ray_origin) / L2_norm(ray_direction)
+
+ # return hit data
+ cdef RayHitData hd
+ hd.u = u
+ hd.v = v
+ hd.t = t
+ hd.converged = (iterations < max_iter)
+ return hd
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_tet_patch_intersect(const void* primitives,
+ const np.int64_t item,
+ Ray* ray) nogil:
+
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ cdef RayHitData hd = compute_tet_patch_hit(tet_patch.v, ray.origin, ray.direction)
+
+ # only count this is it's the closest hit
+ if (hd.t < ray.t_near or hd.t > ray.t_far):
+ return False
+
+ if (0 <= hd.u and 0 <= hd.v and hd.u + hd.v <= 1.0 and hd.converged):
+ # we have a hit, so update ray information
+ ray.t_far = hd.t
+ ray.elem_id = tet_patch.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patch_centroid(const void *primitives,
+ const np.int64_t item,
+ np.float64_t[3] centroid) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ centroid[j] = 0.0
+
+ for i in range(6):
+ for j in range(3):
+ centroid[j] += tet_patch.v[i][j]
+
+ for j in range(3):
+ centroid[j] /= 6.0
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void tet_patch_bbox(const void *primitives,
+ const np.int64_t item,
+ BBox* bbox) nogil:
+
+ cdef np.int64_t i, j
+ cdef TetPatch tet_patch = (<TetPatch*> primitives)[item]
+
+ for j in range(3):
+ bbox.left_edge[j] = tet_patch.v[0][j]
+ bbox.right_edge[j] = tet_patch.v[0][j]
+
+ for i in range(1, 6):
+ for j in range(3):
+ bbox.left_edge[j] = fmin(bbox.left_edge[j], tet_patch.v[i][j])
+ bbox.right_edge[j] = fmax(bbox.right_edge[j], tet_patch.v[i][j])
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/utilities/orientation.py
--- a/yt/utilities/orientation.py
+++ b/yt/utilities/orientation.py
@@ -23,7 +23,7 @@
def _aligned(a, b):
aligned_component = np.abs(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
return np.isclose(aligned_component, 1.0, 1.0e-13)
-
+
def _validate_unit_vectors(normal_vector, north_vector):
@@ -35,7 +35,6 @@
if not np.dot(normal_vector, normal_vector) > 0:
raise YTException("normal_vector cannot be the zero vector.")
-
if north_vector is not None and _aligned(north_vector, normal_vector):
raise YTException("normal_vector and north_vector cannot be aligned.")
@@ -85,7 +84,7 @@
t = np.cross(normal_vector, vecs).sum(axis=1)
ax = t.argmax()
east_vector = np.cross(vecs[ax, :], normal_vector).ravel()
- # self.north_vector must remain None otherwise rotations about a fixed axis will break.
+ # self.north_vector must remain None otherwise rotations about a fixed axis will break.
# The north_vector calculated here will still be included in self.unit_vectors.
north_vector = np.cross(normal_vector, east_vector).ravel()
else:
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -158,7 +158,7 @@
def position():
doc = '''
- The location of the camera.
+ The location of the camera.
Parameters
----------
@@ -179,7 +179,7 @@
'Cannot set the camera focus and position to the same value')
self._position = position
self.switch_orientation(normal_vector=self.focus - self._position,
- north_vector=None)
+ north_vector=self.north_vector)
def fdel(self):
del self._position
@@ -386,9 +386,9 @@
calculated automatically.
"""
+ if north_vector is not None:
+ self.north_vector = north_vector
self.position = position
- self.switch_orientation(normal_vector=self.focus - self.position,
- north_vector=north_vector)
def get_position(self):
"""Return the current camera position"""
@@ -483,11 +483,11 @@
>>> sc = Scene()
>>> cam = sc.add_camera()
>>> # rotate the camera by pi / 4 radians:
- >>> cam.rotate(np.pi/4.0)
+ >>> cam.rotate(np.pi/4.0)
>>> # rotate the camera about the y-axis instead of cam.north_vector:
- >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, np.array([0.0, 1.0, 0.0]))
>>> # rotate the camera about the origin instead of its own position:
- >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
+ >>> cam.rotate(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
"""
rotate_all = rot_vector is not None
@@ -593,7 +593,7 @@
>>> sc = Scene()
>>> cam = sc.add_camera(ds)
>>> # roll the camera by pi / 4 radians:
- >>> cam.roll(np.pi/4.0)
+ >>> cam.roll(np.pi/4.0)
>>> # roll the camera about the origin instead of its own position:
>>> cam.roll(np.pi/4.0, rot_center=np.array([0.0, 0.0, 0.0]))
@@ -627,7 +627,7 @@
>>> import yt
>>> import numpy as np
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
>>> cam = sc.camera
>>> for i in cam.iter_rotate(np.pi, 10):
@@ -688,7 +688,7 @@
def zoom(self, factor):
r"""Change the width of the FOV of the camera.
- This will appear to zoom the camera in by some `factor` toward the
+ This will appear to zoom the camera in by some `factor` toward the
focal point along the current view direction, but really it's just
changing the width of the field of view.
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/render_source.py
--- a/yt/visualization/volume_rendering/render_source.py
+++ b/yt/visualization/volume_rendering/render_source.py
@@ -131,7 +131,7 @@
--------
The easiest way to make a VolumeSource is to use the volume_render
- function, so that the VolumeSource gets created automatically. This
+ function, so that the VolumeSource gets created automatically. This
example shows how to do this and then access the resulting source:
>>> import yt
@@ -545,7 +545,7 @@
'''
This is the name of the colormap that will be used when rendering
this MeshSource object. Should be a string, like 'arbre', or 'dusk'.
-
+
'''
def fget(self):
@@ -562,7 +562,7 @@
'''
These are the bounds that will be used with the colormap to the display
the rendered image. Should be a (vmin, vmax) tuple, like (0.0, 2.0). If
- None, the bounds will be automatically inferred from the max and min of
+ None, the bounds will be automatically inferred from the max and min of
the rendered data.
'''
@@ -603,10 +603,10 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry. Right now, high-order geometry is only
# implemented for 20-point hexes.
- if indices.shape[1] == 20:
+ if indices.shape[1] == 20 or indices.shape[1] == 10:
self.mesh = mesh_construction.QuadraticElementMesh(self.volume,
vertices,
indices,
@@ -620,12 +620,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
self.mesh = mesh_construction.LinearElementMesh(self.volume,
vertices,
@@ -651,7 +645,7 @@
if len(field_data.shape) == 1:
field_data = np.expand_dims(field_data, 1)
- # Here, we decide whether to render based on high-order or
+ # Here, we decide whether to render based on high-order or
# low-order geometry.
if indices.shape[1] == 27:
# hexahedral
@@ -659,12 +653,6 @@
"dropping to 1st order.")
field_data = field_data[:, 0:8]
indices = indices[:, 0:8]
- elif indices.shape[1] == 10:
- # tetrahedral
- mylog.warning("10-node tetrahedral elements not yet supported, " +
- "dropping to 1st order.")
- field_data = field_data[:, 0:4]
- indices = indices[:, 0:4]
self.volume = BVH(vertices, indices, field_data)
@@ -794,7 +782,7 @@
cmap_name=self._cmap)/255.
alpha = image[:, :, 3]
alpha[self.sampler.aimage_used == -1] = 0.0
- image[:, :, 3] = alpha
+ image[:, :, 3] = alpha
return image
def __repr__(self):
@@ -854,7 +842,7 @@
assert(positions.ndim == 2 and positions.shape[1] == 3)
if colors is not None:
assert(colors.ndim == 2 and colors.shape[1] == 4)
- assert(colors.shape[0] == positions.shape[0])
+ assert(colors.shape[0] == positions.shape[0])
self.positions = positions
# If colors aren't individually set, make black with full opacity
if colors is None:
@@ -1057,7 +1045,7 @@
Examples
--------
- This example shows how to use BoxSource to add an outline of the
+ This example shows how to use BoxSource to add an outline of the
domain boundaries to a volume rendering.
>>> import yt
@@ -1065,12 +1053,12 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> box_source = BoxSource(ds.domain_left_edge,
... ds.domain_right_edge,
... [1.0, 1.0, 1.0, 1.0])
>>> sc.add_source(box_source)
- >>>
+ >>>
>>> im = sc.render()
"""
@@ -1078,7 +1066,7 @@
assert(left_edge.shape == (3,))
assert(right_edge.shape == (3,))
-
+
if color is None:
color = np.array([1.0, 1.0, 1.0, 1.0])
@@ -1118,7 +1106,7 @@
Examples
--------
- This example makes a volume rendering and adds outlines of all the
+ This example makes a volume rendering and adds outlines of all the
AMR grids in the simulation:
>>> import yt
@@ -1140,12 +1128,12 @@
>>> import yt
>>> from yt.visualization.volume_rendering.api import GridSource
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
- >>>
+ >>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> dd = ds.sphere("c", (0.1, "unitary"))
>>> grid_source = GridSource(dd, alpha=1.0)
- >>>
+ >>>
>>> sc.add_source(grid_source)
>>>
>>> im = sc.render()
@@ -1223,11 +1211,11 @@
>>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
>>>
>>> im, sc = yt.volume_render(ds)
- >>>
+ >>>
>>> coord_source = CoordinateVectorSource()
- >>>
+ >>>
>>> sc.add_source(coord_source)
- >>>
+ >>>
>>> im = sc.render()
"""
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/tests/test_camera_attributes.py
--- a/yt/visualization/volume_rendering/tests/test_camera_attributes.py
+++ b/yt/visualization/volume_rendering/tests/test_camera_attributes.py
@@ -110,3 +110,12 @@
for lens_type in valid_lens_types:
cam.set_lens(lens_type)
+
+ # See issue #1287
+ cam.focus = [0, 0, 0]
+ cam_pos = [1, 0, 0]
+ north_vector = [0, 1, 0]
+ cam.set_position(cam_pos, north_vector)
+ cam_pos = [0, 1, 0]
+ north_vector = [0, 0, 1]
+ cam.set_position(cam_pos, north_vector)
diff -r f1929d10804e6b316a34b74293ca3c7bdebc26b4 -r 1d43d35f96cf7d7ed8aa8e52222a5a415c4e6215 yt/visualization/volume_rendering/tests/test_mesh_render.py
--- a/yt/visualization/volume_rendering/tests/test_mesh_render.py
+++ b/yt/visualization/volume_rendering/tests/test_mesh_render.py
@@ -55,7 +55,7 @@
images.append(im)
return images
-
+
@requires_module("pyembree")
def test_surface_mesh_render_pyembree():
ytcfg["yt", "ray_tracing_engine"] = "embree"
@@ -145,7 +145,7 @@
im = sc.render()
return compare(ds, im, "%s_render_answers_wedge6_%s_%s" %
(engine, field[0], field[1]))
-
+
@requires_ds(wedge6)
@requires_module("pyembree")
def test_wedge6_render_pyembree():
@@ -157,6 +157,28 @@
for field in wedge6_fields:
yield wedge6_render("yt", field)
+tet10 = "SecondOrderTets/tet10_unstructured_out.e"
+tet10_fields = [('connect1', 'uz')]
+
+def tet10_render(engine, field):
+ ytcfg["yt", "ray_tracing_engine"] = engine
+ ds = data_dir_load(tet10, kwargs={'step':-1})
+ sc = create_scene(ds, field)
+ im = sc.render()
+ return compare(ds, im, "%s_render_answers_tet10_%s_%s" %
+ (engine, field[0], field[1]))
+
+ at requires_ds(tet10)
+ at requires_module("pyembree")
+def test_tet10_render_pyembree():
+ for field in tet10_fields:
+ yield tet10_render("embree", field)
+
+ at requires_ds(tet10)
+def test_tet10_render():
+ for field in tet10_fields:
+ yield tet10_render("yt", field)
+
def perspective_mesh_render(engine):
ytcfg["yt", "ray_tracing_engine"] = engine
ds = data_dir_load(hex8)
@@ -169,7 +191,7 @@
cam.resolution = (800, 800)
im = sc.render()
return compare(ds, im, "%s_perspective_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_perspective_mesh_render_pyembree():
@@ -195,7 +217,7 @@
sc.add_source(ms2)
im = sc.render()
return compare(ds, im, "%s_composite_mesh_render" % engine)
-
+
@requires_ds(hex8)
@requires_module("pyembree")
def test_composite_mesh_render_pyembree():
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list