1# distutils: include_dirs = LIB_DIR 2# distutils: libraries = STD_LIBS 3""" 4Geometry selection routines. 5 6 7 8 9""" 10 11 12import numpy as np 13 14cimport cython 15cimport numpy as np 16cimport oct_visitors 17from cython cimport floating 18from libc.math cimport sqrt 19from libc.stdlib cimport free, malloc 20 21from yt.utilities.lib.bitarray cimport ba_get_value, ba_set_value 22from yt.utilities.lib.fnv_hash cimport c_fnv_hash as fnv_hash 23from yt.utilities.lib.fp_utils cimport fclip, fmax, fmin, iclip, imax, imin 24from yt.utilities.lib.geometry_utils cimport ( 25 bounded_morton_dds, 26 decode_morton_64bit, 27 encode_morton_64bit, 28 morton_neighbors_coarse, 29 morton_neighbors_refined, 30) 31from yt.utilities.lib.grid_traversal cimport sampler_function, walk_volume 32from yt.utilities.lib.volume_container cimport VolumeContainer 33 34from .oct_container cimport Oct, OctreeContainer 35from .oct_visitors cimport cind 36 37 38cdef extern from "math.h": 39 double exp(double x) nogil 40 float expf(float x) nogil 41 long double expl(long double x) nogil 42 double floor(double x) nogil 43 double ceil(double x) nogil 44 double fmod(double x, double y) nogil 45 double log2(double x) nogil 46 long int lrint(double x) nogil 47 double fabs(double x) nogil 48 49# use this as an epsilon test for grids aligned with selector 50# define here to avoid the gil later 51cdef np.float64_t grid_eps = np.finfo(np.float64).eps 52grid_eps = 0.0 53 54cdef inline np.float64_t dot(np.float64_t* v1, 55 np.float64_t* v2) nogil: 56 return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] 57 58cdef inline np.float64_t norm(np.float64_t* v) nogil: 59 return sqrt(dot(v, v)) 60 61# These routines are separated into a couple different categories: 62# 63# * Routines for identifying intersections of an object with a bounding box 64# * Routines for identifying cells/points inside a bounding box that 65# intersect with an object 66# * Routines that speed up some type of geometric calculation 67 68# First, bounding box / object intersection routines. 69# These all respect the interface "dobj" and a set of left_edges, right_edges, 70# sometimes also accepting level and mask information. 71 72@cython.boundscheck(False) 73@cython.wraparound(False) 74@cython.cdivision(True) 75def convert_mask_to_indices(np.ndarray[np.uint8_t, ndim=3, cast=True] mask, 76 int count, int transpose = 0): 77 cdef int i, j, k, cpos 78 cdef np.ndarray[np.int64_t, ndim=2] indices 79 indices = np.zeros((count, 3), dtype='int64') 80 cpos = 0 81 for i in range(mask.shape[0]): 82 for j in range(mask.shape[1]): 83 for k in range(mask.shape[2]): 84 if mask[i, j, k] == 1: 85 if transpose == 1: 86 indices[cpos, 0] = k 87 indices[cpos, 1] = j 88 indices[cpos, 2] = i 89 else: 90 indices[cpos, 0] = i 91 indices[cpos, 1] = j 92 indices[cpos, 2] = k 93 cpos += 1 94 return indices 95 96 97@cython.boundscheck(False) 98@cython.wraparound(False) 99@cython.cdivision(True) 100cdef _mask_fill(np.ndarray[np.float64_t, ndim=1] out, 101 np.int64_t offset, 102 np.ndarray[np.uint8_t, ndim=3, cast=True] mask, 103 np.ndarray[floating, ndim=3] vals): 104 cdef np.int64_t count = 0 105 cdef int i, j, k 106 for i in range(mask.shape[0]): 107 for j in range(mask.shape[1]): 108 for k in range(mask.shape[2]): 109 if mask[i, j, k] == 1: 110 out[offset + count] = vals[i, j, k] 111 count += 1 112 return count 113 114def mask_fill(np.ndarray[np.float64_t, ndim=1] out, 115 np.int64_t offset, 116 np.ndarray[np.uint8_t, ndim=3, cast=True] mask, 117 np.ndarray vals): 118 if vals.dtype == np.float32: 119 return _mask_fill[np.float32_t](out, offset, mask, vals) 120 elif vals.dtype == np.float64: 121 return _mask_fill[np.float64_t](out, offset, mask, vals) 122 else: 123 raise RuntimeError 124 125@cython.cdivision(True) 126@cython.boundscheck(False) 127@cython.wraparound(False) 128def points_in_cells( 129 np.float64_t[:] cx, 130 np.float64_t[:] cy, 131 np.float64_t[:] cz, 132 np.float64_t[:] dx, 133 np.float64_t[:] dy, 134 np.float64_t[:] dz, 135 np.float64_t[:] px, 136 np.float64_t[:] py, 137 np.float64_t[:] pz): 138 # Take a list of cells and particles and calculate which particles 139 # are enclosed within one of the cells. This is used for querying 140 # particle fields on clump/contour objects. 141 # We use brute force since the cells are a relatively unordered collection. 142 143 cdef int p, c, n_p, n_c 144 cdef np.ndarray[np.uint8_t, ndim=1, cast=True] mask 145 146 n_p = px.size 147 n_c = cx.size 148 mask = np.zeros(n_p, dtype="bool") 149 150 for p in range(n_p): 151 for c in range(n_c): 152 if (fabs(px[p] - cx[c]) <= 0.5 * dx[c] and 153 fabs(py[p] - cy[c]) <= 0.5 * dy[c] and 154 fabs(pz[p] - cz[c]) <= 0.5 * dz[c]): 155 mask[p] = True 156 break 157 158 return mask 159 160include "_selection_routines/selector_object.pxi" 161include "_selection_routines/point_selector.pxi" 162include "_selection_routines/sphere_selector.pxi" 163include "_selection_routines/region_selector.pxi" 164include "_selection_routines/cut_region_selector.pxi" 165include "_selection_routines/disk_selector.pxi" 166include "_selection_routines/cutting_plane_selector.pxi" 167include "_selection_routines/slice_selector.pxi" 168include "_selection_routines/ortho_ray_selector.pxi" 169include "_selection_routines/ray_selector.pxi" 170include "_selection_routines/data_collection_selector.pxi" 171include "_selection_routines/ellipsoid_selector.pxi" 172include "_selection_routines/grid_selector.pxi" 173include "_selection_routines/octree_subset_selector.pxi" 174include "_selection_routines/indexed_octree_subset_selector.pxi" 175include "_selection_routines/always_selector.pxi" 176include "_selection_routines/compose_selector.pxi" 177include "_selection_routines/halo_particles_selector.pxi" 178include "_selection_routines/boolean_selectors.pxi" 179