[yt-svn] commit/yt: 6 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jun 8 08:21:28 PDT 2016
6 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/82e8eb2e3b75/
Changeset: 82e8eb2e3b75
Branch: yt
User: chummels
Date: 2016-05-24 23:36:26+00:00
Summary: Adding ion_fraction fields to HI and HII in gizmo frontend
Affected #: 1 file
diff -r 1a33d6bc785673d4cd8ed19beeb2c4de2018a9ae -r 82e8eb2e3b75b0aa90de0f4aea89ea84107c260e yt/frontends/gizmo/fields.py
--- a/yt/frontends/gizmo/fields.py
+++ b/yt/frontends/gizmo/fields.py
@@ -93,6 +93,25 @@
units=self.ds.unit_system["density"])
add_species_field_by_density(self, ptype, "H_p1", particle_type=True)
+ def _h_p0_ion_fraction(field, data):
+ return data[(ptype, 'NeutralHydrogenAbundance')]
+
+ def _h_p1_ion_fraction(field, data):
+ return 1.0 - data[(ptype, 'NeutralHydrogenAbundance')]
+
+ self.add_field(
+ (ptype, "H_ion_fraction"),
+ function=_h_p0_ion_fraction,
+ particle_type=True,
+ units="")
+
+ self.add_field(
+ (ptype, "H_p1_ion_fraction"),
+ function=_h_p1_ion_fraction,
+ particle_type=True,
+ units="")
+ self.alias((ptype, "H_p0_ion_fraction"), (ptype, "H_ion_fraction"))
+
def _nuclei_mass_density_field(field, data):
species = field.name[1][:field.name[1].find("_")]
return data[ptype, "density"] * \
https://bitbucket.org/yt_analysis/yt/commits/cb63dcd8a3fd/
Changeset: cb63dcd8a3fd
Branch: yt
User: chummels
Date: 2016-05-25 01:06:16+00:00
Summary: Assuring hydrogen ion fields get deposited to the grid upon creation.
Affected #: 1 file
diff -r 82e8eb2e3b75b0aa90de0f4aea89ea84107c260e -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 yt/frontends/gizmo/fields.py
--- a/yt/frontends/gizmo/fields.py
+++ b/yt/frontends/gizmo/fields.py
@@ -118,6 +118,15 @@
data[ptype, "%s_metallicity" % species]
num_neighbors = 64
+ for species in ['H', 'H_p0', 'H_p1']:
+ for suf in ["_density", "_mass", "_number_density"]:
+ field = "%s%s" % (species, suf)
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", field,
+ self, num_neighbors)
+ self.alias(("gas", field), fn[0])
+
for species in self.nuclei_names:
self.add_field(
(ptype, "%s_nuclei_mass_density" % species),
https://bitbucket.org/yt_analysis/yt/commits/19a3f7b9bb0b/
Changeset: 19a3f7b9bb0b
Branch: yt
User: chummels
Date: 2016-05-25 06:06:09+00:00
Summary: Merging.
Affected #: 18 files
diff -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 yt/fields/field_info_container.py
--- a/yt/fields/field_info_container.py
+++ b/yt/fields/field_info_container.py
@@ -253,12 +253,17 @@
self[name] = DerivedField(name, f, **kwargs)
return f
return create_function
- ptype = kwargs.get("particle_type", False)
- if ptype:
+
+ if isinstance(name, tuple):
+ self[name] = DerivedField(name, function, **kwargs)
+ return
+
+ if kwargs.get("particle_type", False):
ftype = 'all'
else:
ftype = self.ds.default_fluid_type
- if not isinstance(name, tuple) and (ftype, name) not in self:
+
+ if (ftype, name) not in self:
tuple_name = (ftype, name)
self[tuple_name] = DerivedField(tuple_name, function, **kwargs)
self.alias(name, tuple_name)
diff -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -57,10 +57,6 @@
mylog.debug("Detecting fields.")
self._detect_output_fields()
- def __del__(self):
- if self._data_file is not None:
- self._data_file.close()
-
def _initialize_state_variables(self):
self._parallel_locking = False
self._data_file = None
diff -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 yt/utilities/file_handler.py
--- a/yt/utilities/file_handler.py
+++ b/yt/utilities/file_handler.py
@@ -14,7 +14,6 @@
#-----------------------------------------------------------------------------
from yt.utilities.on_demand_imports import _h5py as h5py
-from distutils.version import LooseVersion
class HDF5FileHandler(object):
handle = None
@@ -22,15 +21,6 @@
def __init__(self, filename):
self.handle = h5py.File(filename, 'r')
- def __del__(self):
- if self.handle is not None:
- # In h5py 2.4 and newer file handles are closed automatically.
- # so only close the handle on old versions. This works around an
- # issue in h5py when run under python 3.4.
- # See https://github.com/h5py/h5py/issues/534
- if LooseVersion(h5py.__version__) < '2.4.0':
- self.handle.close()
-
def __getitem__(self, key):
return self.handle[key]
diff -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 yt/utilities/on_demand_imports.py
--- a/yt/utilities/on_demand_imports.py
+++ b/yt/utilities/on_demand_imports.py
@@ -10,6 +10,8 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+from pkg_resources import parse_version
+
class NotAModule(object):
"""
A class to implement an informative error message that will be outputted if
@@ -28,7 +30,7 @@
def __call__(self, *args, **kwargs):
raise self.error
-class astropy_imports:
+class astropy_imports(object):
_name = "astropy"
_pyfits = None
@property
@@ -105,7 +107,7 @@
_astropy = astropy_imports()
-class scipy_imports:
+class scipy_imports(object):
_name = "scipy"
_integrate = None
@property
@@ -186,12 +188,27 @@
_scipy = scipy_imports()
-class h5py_imports:
+class h5py_imports(object):
_name = "h5py"
+ _err = None
+
+ def __init__(self):
+ try:
+ import h5py
+ if parse_version(h5py.__version__) < parse_version('2.4.0'):
+ self._err = RuntimeError(
+ 'yt requires h5py version 2.4.0 or newer, '
+ 'please update h5py with e.g. "pip install -U h5py" '
+ 'and try again')
+ except ImportError:
+ pass
+ super(h5py_imports, self).__init__()
_File = None
@property
def File(self):
+ if self._err:
+ raise self._err
if self._File is None:
try:
from h5py import File
@@ -203,6 +220,8 @@
_Group = None
@property
def Group(self):
+ if self._err:
+ raise self._err
if self._Group is None:
try:
from h5py import Group
@@ -214,6 +233,8 @@
___version__ = None
@property
def __version__(self):
+ if self._err:
+ raise self._err
if self.___version__ is None:
try:
from h5py import __version__
@@ -225,6 +246,8 @@
_get_config = None
@property
def get_config(self):
+ if self._err:
+ raise self._err
if self._get_config is None:
try:
from h5py import get_config
@@ -236,6 +259,8 @@
_h5f = None
@property
def h5f(self):
+ if self._err:
+ raise self._err
if self._h5f is None:
try:
import h5py.h5f as h5f
@@ -247,6 +272,8 @@
_h5d = None
@property
def h5d(self):
+ if self._err:
+ raise self._err
if self._h5d is None:
try:
import h5py.h5d as h5d
@@ -258,6 +285,8 @@
_h5s = None
@property
def h5s(self):
+ if self._err:
+ raise self._err
if self._h5s is None:
try:
import h5py.h5s as h5s
@@ -269,6 +298,8 @@
_version = None
@property
def version(self):
+ if self._err:
+ raise self._err
if self._version is None:
try:
import h5py.version as version
@@ -279,7 +310,7 @@
_h5py = h5py_imports()
-class nose_imports:
+class nose_imports(object):
_name = "nose"
_run = None
@property
diff -r cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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 cb63dcd8a3fd44edaae2c4e4bfcb2f5d641cfdc0 -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 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)
+
https://bitbucket.org/yt_analysis/yt/commits/d29d9be3bc63/
Changeset: d29d9be3bc63
Branch: yt
User: chummels
Date: 2016-05-25 06:08:12+00:00
Summary: Removing ion_fraction field definition.
Affected #: 1 file
diff -r 19a3f7b9bb0bd3db4fa7cb1b75273a77ae407718 -r d29d9be3bc63be45ad287a0bf29075dd66d18667 yt/frontends/gizmo/fields.py
--- a/yt/frontends/gizmo/fields.py
+++ b/yt/frontends/gizmo/fields.py
@@ -93,25 +93,6 @@
units=self.ds.unit_system["density"])
add_species_field_by_density(self, ptype, "H_p1", particle_type=True)
- def _h_p0_ion_fraction(field, data):
- return data[(ptype, 'NeutralHydrogenAbundance')]
-
- def _h_p1_ion_fraction(field, data):
- return 1.0 - data[(ptype, 'NeutralHydrogenAbundance')]
-
- self.add_field(
- (ptype, "H_ion_fraction"),
- function=_h_p0_ion_fraction,
- particle_type=True,
- units="")
-
- self.add_field(
- (ptype, "H_p1_ion_fraction"),
- function=_h_p1_ion_fraction,
- particle_type=True,
- units="")
- self.alias((ptype, "H_p0_ion_fraction"), (ptype, "H_ion_fraction"))
-
def _nuclei_mass_density_field(field, data):
species = field.name[1][:field.name[1].find("_")]
return data[ptype, "density"] * \
https://bitbucket.org/yt_analysis/yt/commits/b85c646d1012/
Changeset: b85c646d1012
Branch: yt
User: chummels
Date: 2016-06-01 05:40:47+00:00
Summary: Removing H ion mass from the list of deposited fields in gizmo frontend.
Affected #: 1 file
diff -r d29d9be3bc63be45ad287a0bf29075dd66d18667 -r b85c646d10125cd8362e99f25877e9fb3ab45ae0 yt/frontends/gizmo/fields.py
--- a/yt/frontends/gizmo/fields.py
+++ b/yt/frontends/gizmo/fields.py
@@ -100,7 +100,7 @@
num_neighbors = 64
for species in ['H', 'H_p0', 'H_p1']:
- for suf in ["_density", "_mass", "_number_density"]:
+ for suf in ["_density", "_number_density"]:
field = "%s%s" % (species, suf)
fn = add_volume_weighted_smoothed_field(
ptype, "particle_position", "particle_mass",
https://bitbucket.org/yt_analysis/yt/commits/19438cc7e51d/
Changeset: 19438cc7e51d
Branch: yt
User: ngoldbaum
Date: 2016-06-08 15:21:17+00:00
Summary: Merged in chummels/yt (pull request #2195)
Deposit hydrogen fields to grid in gizmo frontend
Affected #: 1 file
diff -r ce2031d2484da5d02bcced20938161f624fccf6f -r 19438cc7e51d035905e59b765a5031836699534d yt/frontends/gizmo/fields.py
--- a/yt/frontends/gizmo/fields.py
+++ b/yt/frontends/gizmo/fields.py
@@ -99,6 +99,15 @@
data[ptype, "%s_metallicity" % species]
num_neighbors = 64
+ for species in ['H', 'H_p0', 'H_p1']:
+ for suf in ["_density", "_number_density"]:
+ field = "%s%s" % (species, suf)
+ fn = add_volume_weighted_smoothed_field(
+ ptype, "particle_position", "particle_mass",
+ "smoothing_length", "density", field,
+ self, num_neighbors)
+ self.alias(("gas", field), fn[0])
+
for species in self.nuclei_names:
self.add_field(
(ptype, "%s_nuclei_mass_density" % species),
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