[yt-svn] commit/yt: ngoldbaum: Merged in atmyers/yt (pull request #2174)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon May 23 20:30:32 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/19aead8b2225/
Changeset: 19aead8b2225
Branch: yt
User: ngoldbaum
Date: 2016-05-24 03:30:21+00:00
Summary: Merged in atmyers/yt (pull request #2174)
Changing the way mesh line annotations are performed.
Affected #: 14 files
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 doc/source/visualizing/plots.rst
--- a/doc/source/visualizing/plots.rst
+++ b/doc/source/visualizing/plots.rst
@@ -406,16 +406,14 @@
import yt
ds = yt.load('MOOSE_sample_data/out.e-s010')
sl = yt.SlicePlot(ds, 'z', ('connect1', 'diffused'))
- sl.annotate_mesh_lines(thresh=0.1)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
sl.zoom(0.75)
sl.save()
-This annotation is performed by marking the pixels where the mapped coordinate is close
-to the element boundary. What counts as 'close' (in the mapped coordinate system) is
-determined by the ``thresh`` parameter, which can be varied to make the lines thicker or
-thinner.
+The ``plot_args`` parameter is a dictionary of keyword arguments that will be passed
+to matplotlib. It can be used to control the mesh line color, thickness, etc...
-The above example all involve 8-node hexahedral mesh elements. Here is another example from
+The above examples all involve 8-node hexahedral mesh elements. Here is another example from
a dataset that uses 6-node wedge elements:
.. python-script::
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 setup.py
--- a/setup.py
+++ b/setup.py
@@ -157,12 +157,6 @@
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs, depends=["yt/utilities/lib/bitarray.pxd"]),
- Extension("yt.utilities.lib.primitives",
- ["yt/utilities/lib/primitives.pyx"],
- libraries=std_libs,
- depends=["yt/utilities/lib/primitives.pxd",
- "yt/utilities/lib/vec3_ops.pxd",
- "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.bounding_volume_hierarchy",
["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
include_dirs=["yt/utilities/lib/"],
@@ -170,6 +164,7 @@
extra_link_args=omp_args,
libraries=std_libs,
depends=["yt/utilities/lib/element_mappings.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/vec3_ops.pxd",
"yt/utilities/lib/primitives.pxd"]),
Extension("yt.utilities.lib.contour_finding",
@@ -196,14 +191,22 @@
"yt/utilities/lib/fixed_interpolator.pxd",
"yt/utilities/lib/fixed_interpolator.h",
]),
+ Extension("yt.utilities.lib.mesh_triangulation",
+ ["yt/utilities/lib/mesh_triangulation.pyx"],
+ depends=["yt/utilities/lib/mesh_triangulation.h"]),
Extension("yt.utilities.lib.pixelization_routines",
["yt/utilities/lib/pixelization_routines.pyx",
"yt/utilities/lib/pixelization_constants.c"],
include_dirs=["yt/utilities/lib/"],
- language="c++",
libraries=std_libs, depends=["yt/utilities/lib/fp_utils.pxd",
"yt/utilities/lib/pixelization_constants.h",
"yt/utilities/lib/element_mappings.pxd"]),
+ Extension("yt.utilities.lib.primitives",
+ ["yt/utilities/lib/primitives.pyx"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/primitives.pxd",
+ "yt/utilities/lib/vec3_ops.pxd",
+ "yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.origami",
["yt/utilities/lib/origami.pyx",
"yt/utilities/lib/origami_tags.c"],
@@ -283,6 +286,7 @@
Extension("yt.utilities.lib.mesh_construction",
["yt/utilities/lib/mesh_construction.pyx"],
depends=["yt/utilities/lib/mesh_construction.pxd",
+ "yt/utilities/lib/mesh_triangulation.h",
"yt/utilities/lib/mesh_intersection.pxd",
"yt/utlilites/lib/mesh_samplers.pxd",
"yt/utlilites/lib/mesh_traversal.pxd"]),
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -76,10 +76,15 @@
field_data = ad[field]
buff_size = size[0:dimension] + (1,) + size[dimension:]
+ ax = data_source.axis
+ xax = self.x_axis[ax]
+ yax = self.y_axis[ax]
c = np.float64(data_source.center[dimension].d)
- bounds.insert(2*dimension, c)
- bounds.insert(2*dimension, c)
- bounds = np.reshape(bounds, (3, 2))
+
+ extents = np.zeros((3, 2))
+ extents[ax] = np.array([c, c])
+ extents[xax] = bounds[0:2]
+ extents[yax] = bounds[2:4]
# if this is an element field, promote to 2D here
if len(field_data.shape) == 1:
@@ -101,13 +106,10 @@
img = pixelize_element_mesh(coords,
indices,
- buff_size, field_data, bounds,
+ buff_size, field_data, extents,
index_offset=offset)
# re-order the array and squeeze out the dummy dim
- ax = data_source.axis
- xax = self.x_axis[ax]
- yax = self.y_axis[ax]
return np.squeeze(np.transpose(img, (yax, xax, ax)))
elif dimension < 3:
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- a/yt/utilities/lib/bounding_volume_hierarchy.pxd
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -4,7 +4,7 @@
from yt.utilities.lib.element_mappings cimport ElementSampler
from yt.utilities.lib.primitives cimport BBox, Ray
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -30,11 +30,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -52,11 +47,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef class P1Sampler3D(ElementSampler):
@@ -72,11 +62,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -169,11 +154,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -191,11 +171,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -214,11 +189,6 @@
cdef int check_inside(self, double* mapped_coord) nogil
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
-
cdef int check_mesh_lines(self, double* mapped_coord) nogil
@@ -250,8 +220,3 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
-
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -86,15 +86,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- pass
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
pass
@@ -182,19 +173,6 @@
return 0
return 1
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
cdef class P1Sampler3D(ElementSampler):
'''
@@ -261,19 +239,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
-
- if (fabs(mapped_coord[direction]) < tolerance or
- fabs(mapped_coord[direction] - 1.0) < tolerance):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double u, v, w
cdef double thresh = 2.0e-2
@@ -413,18 +378,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
if (fabs(fabs(mapped_coord[0]) - 1.0) < 1e-1 and
fabs(fabs(mapped_coord[1]) - 1.0) < 1e-1):
@@ -568,19 +521,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -781,22 +721,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (direction == 2 and (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance)):
- return 1
- elif (direction == 1 and (fabs(mapped_coord[direction]) < tolerance)):
- return 1
- elif (direction == 0 and (fabs(mapped_coord[0]) < tolerance or
- fabs(mapped_coord[0] + mapped_coord[1] - 1.0) < tolerance)):
- return 1
- return 0
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
cdef int check_mesh_lines(self, double* mapped_coord) nogil:
cdef double r, s, t
cdef double thresh = 5.0e-2
@@ -969,20 +893,6 @@
return 0
return 1
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef int check_near_edge(self,
- double* mapped_coord,
- double tolerance,
- int direction) nogil:
- if (fabs(fabs(mapped_coord[direction]) - 1.0) < tolerance):
- return 1
- else:
- return 0
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/geometry_utils.pyx
--- a/yt/utilities/lib/geometry_utils.pyx
+++ b/yt/utilities/lib/geometry_utils.pyx
@@ -18,8 +18,9 @@
cimport cython
from libc.stdlib cimport malloc, free
from yt.utilities.lib.fp_utils cimport fclip, i64clip
-from libc.math cimport copysign
+from libc.math cimport copysign, fabs
from yt.utilities.exceptions import YTDomainOverflow
+from yt.utilities.lib.vec3_ops cimport subtract, cross, dot, L2_norm
cdef extern from "math.h":
double exp(double x) nogil
@@ -501,6 +502,11 @@
cdef np.float64_t p0[3]
cdef np.float64_t p1[3]
cdef np.float64_t p3[3]
+ cdef np.float64_t E0[3]
+ cdef np.float64_t E1[3]
+ cdef np.float64_t tri_norm[3]
+ cdef np.float64_t plane_norm[3]
+ cdef np.float64_t dp
cdef int i, j, k, count, ntri, nlines
nlines = 0
ntri = triangles.shape[0]
@@ -510,6 +516,22 @@
first = last = points = NULL
for i in range(ntri):
count = 0
+
+ # skip if triangle is close to being parallel to plane
+ for j in range(3):
+ p0[j] = triangles[i, 0, j]
+ p1[j] = triangles[i, 1, j]
+ p3[j] = triangles[i, 2, j]
+ plane_norm[j] = 0.0
+ plane_norm[ax] = 1.0
+ subtract(p1, p0, E0)
+ subtract(p3, p0, E1)
+ cross(E0, E1, tri_norm)
+ dp = dot(tri_norm, plane_norm)
+ dp /= L2_norm(tri_norm)
+ if (fabs(dp) > 0.995):
+ continue
+
# Now for each line segment (01, 12, 20) we check to see how many cross
# the coordinate.
for j in range(3):
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_construction.h
--- a/yt/utilities/lib/mesh_construction.h
+++ /dev/null
@@ -1,66 +0,0 @@
-#define MAX_NUM_TRI 12
-#define HEX_NV 8
-#define HEX_NT 12
-#define TETRA_NV 4
-#define TETRA_NT 4
-#define WEDGE_NV 6
-#define WEDGE_NT 8
-
-// This array is used to triangulate the hexahedral mesh elements
-// Each element has six faces with two triangles each.
-// The vertex ordering convention is assumed to follow that used
-// 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
- {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
- {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
- {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
-};
-
-// Similarly, this is used to triangulate the tetrahedral cells
-int triangulate_tetra[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 1, 3},
- {0, 2, 3},
- {1, 2, 3},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-// Triangulate wedges
-int triangulate_wedge[MAX_NUM_TRI][3] = {
- {0, 1, 2},
- {0, 3, 1},
- {1, 3, 4},
- {0, 2, 3},
- {2, 5, 3},
- {1, 4, 2},
- {2, 4, 5},
- {3, 5, 4},
-
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1},
- {-1, -1, -1}
-};
-
-
-// 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},
- {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}
-};
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_construction.pyx
--- a/yt/utilities/lib/mesh_construction.pyx
+++ b/yt/utilities/lib/mesh_construction.pyx
@@ -40,7 +40,7 @@
patchBoundsFunc
from yt.utilities.exceptions import YTElementTypeNotRecognized
-cdef extern from "mesh_construction.h":
+cdef extern from "mesh_triangulation.h":
enum:
MAX_NUM_TRI
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_triangulation.h
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.h
@@ -0,0 +1,66 @@
+#define MAX_NUM_TRI 12
+#define HEX_NV 8
+#define HEX_NT 12
+#define TETRA_NV 4
+#define TETRA_NT 4
+#define WEDGE_NV 6
+#define WEDGE_NT 8
+
+// This array is used to triangulate the hexahedral mesh elements
+// Each element has six faces with two triangles each.
+// The vertex ordering convention is assumed to follow that used
+// 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
+ {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
+ {0, 7, 3}, {0, 4, 7}, // Face is 3 0 4 7
+ {3, 6, 2}, {3, 7, 6} // Face is 2 3 7 6
+};
+
+// Similarly, this is used to triangulate the tetrahedral cells
+int triangulate_tetra[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 1, 3},
+ {0, 2, 3},
+ {1, 2, 3},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+// Triangulate wedges
+int triangulate_wedge[MAX_NUM_TRI][3] = {
+ {0, 1, 2},
+ {0, 3, 1},
+ {1, 3, 4},
+ {0, 2, 3},
+ {2, 5, 3},
+ {1, 4, 2},
+ {2, 4, 5},
+ {3, 5, 4},
+
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1},
+ {-1, -1, -1}
+};
+
+
+// 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},
+ {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}
+};
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/mesh_triangulation.pyx
--- /dev/null
+++ b/yt/utilities/lib/mesh_triangulation.pyx
@@ -0,0 +1,58 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+
+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
+ 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]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def triangulate_indices(np.ndarray[np.int64_t, ndim=2] indices):
+ cdef np.int64_t num_elem = indices.shape[0]
+ cdef np.int64_t VPE = indices.shape[1] # num verts per element
+ cdef np.int64_t TPE # num triangles per element
+ cdef int[MAX_NUM_TRI][3] tri_array
+
+ if (VPE == 8 or VPE == 20 or VPE == 27):
+ TPE = HEX_NT
+ tri_array = triangulate_hex
+ elif VPE == 4:
+ TPE = TETRA_NT
+ tri_array = triangulate_tetra
+ elif VPE == 6:
+ TPE = WEDGE_NT
+ tri_array = triangulate_wedge
+ else:
+ raise YTElementTypeNotRecognized(3, VPE)
+
+ cdef np.int64_t num_tri = TPE * num_elem
+
+ cdef np.ndarray[np.int64_t, ndim=2] tri_indices
+ tri_indices = np.empty((num_tri, 3), dtype="int64")
+
+ cdef np.int64_t *tri_indices_ptr = <np.int64_t*> tri_indices.data
+ cdef np.int64_t *indices_ptr = <np.int64_t*> indices.data
+
+ cdef np.int64_t i, j, k, offset
+ for i in range(num_elem):
+ for j in range(TPE):
+ for k in range(3):
+ offset = tri_array[j][k]
+ tri_indices_ptr[i*TPE*3 + 3*j + k] = indices_ptr[i*VPE + offset]
+ return tri_indices
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -18,7 +18,8 @@
cimport cython
cimport libc.math as math
from yt.utilities.lib.fp_utils cimport fmin, fmax, i64min, i64max, imin, imax, fabs
-from yt.utilities.exceptions import YTPixelizeError, \
+from yt.utilities.exceptions import \
+ YTPixelizeError, \
YTElementTypeNotRecognized
from libc.stdlib cimport malloc, free
from vec3_ops cimport dot, cross, subtract
@@ -544,8 +545,7 @@
buff_size,
np.ndarray[np.float64_t, ndim=2] field,
extents,
- int index_offset = 0,
- double thresh = -1.0):
+ int index_offset = 0):
cdef np.ndarray[np.float64_t, ndim=3] img
img = np.zeros(buff_size, dtype="float64")
# Two steps:
@@ -592,27 +592,10 @@
# if we are in 2D land, the 1 cell thick dimension had better be 'z'
if ndim == 2:
- assert(buff_size[2] == 1)
+ if buff_size[2] != 1:
+ raise RuntimeError("Slices of 2D datasets must be "
+ "perpendicular to the 'z' direction.")
- ax = -1
- for i in range(3):
- if buff_size[i] == 1:
- ax = i
- if ax == -1:
- raise RuntimeError
- xax = yax = -1
- if ax == 0:
- xax = 1
- yax = 2
- elif ax == 1:
- xax = 1
- yax = 2
- elif ax == 2:
- xax = 0
- yax = 1
- if xax == -1 or yax == -1:
- raise RuntimeError
-
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
@@ -673,20 +656,11 @@
sampler.map_real_to_unit(mapped_coord, vertices, ppoint)
if not sampler.check_inside(mapped_coord):
continue
- # if thresh is negative, we do the normal sample operation
- if thresh < 0.0:
- if (num_field_vals == 1):
- img[pi, pj, pk] = field_vals[0]
- else:
- img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
- field_vals)
+ if (num_field_vals == 1):
+ img[pi, pj, pk] = field_vals[0]
else:
- # otherwise, we draw the element boundaries
- if sampler.check_near_edge(mapped_coord, thresh, xax):
- img[pi, pj, pk] = 1.0
- elif sampler.check_near_edge(mapped_coord, thresh, yax):
- img[pi, pj, pk] = 1.0
-
+ img[pi, pj, pk] = sampler.sample_at_unit_point(mapped_coord,
+ field_vals)
free(vertices)
free(field_vals)
return img
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -22,8 +22,6 @@
from distutils.version import LooseVersion
-from yt.config import \
- ytcfg
from yt.funcs import \
mylog, iterable
from yt.extern.six import add_metaclass
@@ -31,13 +29,14 @@
from yt.visualization.image_writer import apply_colormap
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.utilities.lib.pixelization_routines import \
- pixelize_element_mesh, pixelize_off_axis_cartesian, \
+ pixelize_off_axis_cartesian, \
pixelize_cartesian
from yt.analysis_modules.cosmological_observation.light_ray.light_ray import \
periodic_ray
from yt.utilities.lib.line_integral_convolution import \
line_integral_convolution_2d
-
+from yt.geometry.unstructured_mesh_handler import UnstructuredIndex
+from yt.utilities.lib.mesh_triangulation import triangulate_indices
from . import _MPL
@@ -1605,85 +1604,77 @@
class MeshLinesCallback(PlotCallback):
"""
- annotate_mesh_lines(thresh=0.01)
+ annotate_mesh_lines()
- Adds the mesh lines to the plot.
-
- This is done by marking the pixels where the mapped coordinate
- is within some threshold distance of one of the element boundaries.
- If the mesh lines are too thick or too thin, try varying thresh.
+ Adds mesh lines to the plot. Only works for unstructured or
+ semi-structured mesh data. For structured grid data, see
+ GridBoundaryCallback or CellEdgesCallback.
Parameters
----------
- thresh : float
- The threshold distance, in mapped coordinates, within which the
- pixels will be marked as part of the element boundary.
- Default is 0.01.
+ plot_args: dict, optional
+ A dictionary of arguments that will be passed to matplotlib.
+
+ Example
+ -------
+
+ >>> import yt
+ >>> ds = yt.load("MOOSE_sample_data/out.e-s010")
+ >>> sl = yt.SlicePlot(ds, 'z', ('connect2', 'convected'))
+ >>> sl.annotate_mesh_lines(plot_args={'color':'black'})
"""
_type_name = "mesh_lines"
- def __init__(self, thresh=0.1, cmap=None):
+ def __init__(self, plot_args=None):
super(MeshLinesCallback, self).__init__()
- self.thresh = thresh
- if cmap is None:
- cmap = ytcfg.get("yt", "default_colormap")
- self.cmap = cmap
+ self.plot_args = plot_args
+
+ def promote_2d_to_3d(self, coords, indices, plot):
+ new_coords = np.zeros((2*coords.shape[0], 3))
+ new_connects = np.zeros((indices.shape[0], 2*indices.shape[1]),
+ dtype=np.int64)
+
+ new_coords[0:coords.shape[0],0:2] = coords
+ new_coords[0:coords.shape[0],2] = plot.ds.domain_left_edge[2]
+ new_coords[coords.shape[0]:,0:2] = coords
+ new_coords[coords.shape[0]:,2] = plot.ds.domain_right_edge[2]
+
+ new_connects[:,0:indices.shape[1]] = indices
+ new_connects[:,indices.shape[1]:] = indices + coords.shape[0]
+
+ return new_coords, new_connects
def __call__(self, plot):
- ftype, fname = plot.field
- mesh_id = int(ftype[-1]) - 1
- mesh = plot.ds.index.meshes[mesh_id]
+ index = plot.ds.index
+ if not issubclass(type(index), UnstructuredIndex):
+ raise RuntimeError("Mesh line annotations only work for "
+ "unstructured or semi-structured mesh data.")
+ for i, m in enumerate(index.meshes):
+ try:
+ ftype, fname = plot.field
+ if ftype.startswith('connect') and int(ftype[-1]) - 1 != i:
+ continue
+ except ValueError:
+ pass
+ coords = m.connectivity_coords
+ indices = m.connectivity_indices - m._index_offset
+
+ num_verts = indices.shape[1]
+ num_dims = coords.shape[1]
- coords = mesh.connectivity_coords
- indices = mesh.connectivity_indices
+ if num_dims == 2 and num_verts == 3:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
+ elif num_dims == 2 and num_verts == 4:
+ coords, indices = self.promote_2d_to_3d(coords, indices, plot)
- xx0, xx1 = plot._axes.get_xlim()
- yy0, yy1 = plot._axes.get_ylim()
-
- offset = mesh._index_offset
- ax = plot.data.axis
- xax = plot.data.ds.coordinates.x_axis[ax]
- yax = plot.data.ds.coordinates.y_axis[ax]
-
- size = plot.image._A.shape
- c = plot.data.center[ax]
- buff_size = size[0:ax] + (1,) + size[ax:]
-
- x0, x1 = plot.xlim
- y0, y1 = plot.ylim
- c = plot.data.center[ax]
- bounds = [x0, x1, y0, y1]
- bounds.insert(2*ax, c)
- bounds.insert(2*ax, c)
- bounds = np.reshape(bounds, (3, 2))
-
- # draw the mesh lines by marking where the mapped
- # coordinate is close to the domain edge. We do this by
- # calling the pixelizer with a dummy field and passing in
- # a non-negative thresh.
- dummy_field = np.zeros(indices.shape, dtype=np.float64)
- img = pixelize_element_mesh(coords, indices,
- buff_size,
- dummy_field,
- bounds,
- offset, self.thresh)
- img = np.squeeze(np.transpose(img, (yax, xax, ax)))
-
- # convert to RGBA
- size = plot.frb.buff_size
- image = np.zeros((size[0], size[1], 4), dtype=np.uint8)
- image[:, :, 0][img > 0.0] = 0
- image[:, :, 1][img > 0.0] = 0
- image[:, :, 2][img > 0.0] = 0
- image[:, :, 3][img > 0.0] = 255
-
- plot._axes.imshow(image, zorder=1,
- extent=[xx0, xx1, yy0, yy1],
- origin='lower', cmap=self.cmap,
- interpolation='nearest')
+ tri_indices = triangulate_indices(indices)
+ points = coords[tri_indices]
+
+ tfc = TriangleFacetsCallback(points, plot_args=self.plot_args)
+ tfc(plot)
class TriangleFacetsCallback(PlotCallback):
diff -r 8bd28805baeed4b0f741e06810c156717c838579 -r 19aead8b2225bdd19398c18c07470d084f8fbca0 yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -21,7 +21,9 @@
from yt.config import \
ytcfg
from yt.testing import \
- fake_amr_ds
+ fake_amr_ds, \
+ fake_tetrahedral_ds, \
+ fake_hexahedral_ds
import yt.units as u
from .test_plotwindow import assert_fname
from yt.utilities.exceptions import \
@@ -368,6 +370,22 @@
color=(0.0, 1.0, 1.0))
p.save(prefix)
+def test_mesh_lines_callback():
+ with _cleanup_fname() as prefix:
+
+ ds = fake_hexahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+ ds = fake_tetrahedral_ds()
+ for field in ds.field_list:
+ sl = SlicePlot(ds, 1, field)
+ sl.annotate_mesh_lines(plot_args={'color':'black'})
+ yield assert_fname, sl.save(prefix)[0]
+
+
def test_line_integral_convolution_callback():
with _cleanup_fname() as prefix:
ds = fake_amr_ds(fields =
@@ -390,3 +408,4 @@
alpha=0.9, const_alpha=True)
p.save(prefix)
+
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