1cdef class SphereSelector(SelectorObject): 2 cdef public np.float64_t radius 3 cdef public np.float64_t radius2 4 cdef public np.float64_t center[3] 5 cdef np.float64_t bbox[3][2] 6 cdef public bint check_box[3] 7 8 def __init__(self, dobj): 9 for i in range(3): 10 self.center[i] = _ensure_code(dobj.center[i]) 11 self.radius = _ensure_code(dobj.radius) 12 self.radius2 = self.radius * self.radius 13 self.set_bbox(_ensure_code(dobj.center)) 14 for i in range(3): 15 if self.bbox[i][0] < dobj.ds.domain_left_edge[i]: 16 self.check_box[i] = False 17 elif self.bbox[i][1] > dobj.ds.domain_right_edge[i]: 18 self.check_box[i] = False 19 else: 20 self.check_box[i] = True 21 22 def set_bbox(self, center): 23 for i in range(3): 24 self.center[i] = center[i] 25 self.bbox[i][0] = self.center[i] - self.radius 26 self.bbox[i][1] = self.center[i] + self.radius 27 28 @cython.boundscheck(False) 29 @cython.wraparound(False) 30 @cython.cdivision(True) 31 cdef int select_cell(self, np.float64_t pos[3], np.float64_t dds[3]) nogil: 32 # sphere center either inside cell or center of cell lies inside sphere 33 if (pos[0] - 0.5*dds[0] <= self.center[0] <= pos[0]+0.5*dds[0] and 34 pos[1] - 0.5*dds[1] <= self.center[1] <= pos[1]+0.5*dds[1] and 35 pos[2] - 0.5*dds[2] <= self.center[2] <= pos[2]+0.5*dds[2]): 36 return 1 37 return self.select_point(pos) 38 # # langmm: added to allow sphere to interesect edge/corner of cell 39 # cdef np.float64_t LE[3] 40 # cdef np.float64_t RE[3] 41 # cdef int i 42 # for i in range(3): 43 # LE[i] = pos[i] - 0.5*dds[i] 44 # RE[i] = pos[i] + 0.5*dds[i] 45 # return self.select_bbox(LE, RE) 46 47 @cython.boundscheck(False) 48 @cython.wraparound(False) 49 @cython.cdivision(True) 50 cdef int select_point(self, np.float64_t pos[3]) nogil: 51 cdef int i 52 cdef np.float64_t dist, dist2 = 0 53 for i in range(3): 54 if self.check_box[i] and \ 55 (pos[i] < self.bbox[i][0] or 56 pos[i] > self.bbox[i][1]): 57 return 0 58 dist = _periodic_dist(pos[i], self.center[i], self.domain_width[i], 59 self.periodicity[i]) 60 dist2 += dist*dist 61 if dist2 > self.radius2: return 0 62 return 1 63 64 @cython.boundscheck(False) 65 @cython.wraparound(False) 66 @cython.cdivision(True) 67 cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil: 68 cdef int i 69 cdef np.float64_t dist, dist2 = 0 70 for i in range(3): 71 dist = self.periodic_difference(pos[i], self.center[i], i) 72 dist2 += dist*dist 73 dist = self.radius+radius 74 if dist2 <= dist*dist: return 1 75 return 0 76 77 @cython.boundscheck(False) 78 @cython.wraparound(False) 79 @cython.cdivision(True) 80 cdef int select_bbox(self, np.float64_t left_edge[3], 81 np.float64_t right_edge[3]) nogil: 82 cdef np.float64_t box_center, relcenter, closest, dist, edge 83 cdef int i 84 if (left_edge[0] <= self.center[0] < right_edge[0] and 85 left_edge[1] <= self.center[1] < right_edge[1] and 86 left_edge[2] <= self.center[2] < right_edge[2]): 87 return 1 88 for i in range(3): 89 if not self.check_box[i]: continue 90 if right_edge[i] < self.bbox[i][0] or \ 91 left_edge[i] > self.bbox[i][1]: 92 return 0 93 # http://www.gamedev.net/topic/335465-is-this-the-simplest-sphere-aabb-collision-test/ 94 dist = 0 95 for i in range(3): 96 # Early terminate 97 box_center = (right_edge[i] + left_edge[i])/2.0 98 relcenter = self.periodic_difference(box_center, self.center[i], i) 99 edge = right_edge[i] - left_edge[i] 100 closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0) 101 dist += closest*closest 102 if dist > self.radius2: return 0 103 return 1 104 105 @cython.boundscheck(False) 106 @cython.wraparound(False) 107 @cython.cdivision(True) 108 cdef int select_bbox_edge(self, np.float64_t left_edge[3], 109 np.float64_t right_edge[3]) nogil: 110 cdef np.float64_t box_center, relcenter, closest, farthest, cdist, fdist, edge 111 cdef int i 112 if (left_edge[0] <= self.center[0] <= right_edge[0] and 113 left_edge[1] <= self.center[1] <= right_edge[1] and 114 left_edge[2] <= self.center[2] <= right_edge[2]): 115 fdist = 0 116 for i in range(3): 117 edge = right_edge[i] - left_edge[i] 118 box_center = (right_edge[i] + left_edge[i])/2.0 119 relcenter = self.periodic_difference( 120 box_center, self.center[i], i) 121 if relcenter >= 0: 122 farthest = relcenter + edge/2.0 123 else: 124 farthest = relcenter - edge/2.0 125 # farthest = relcenter + fclip(relcenter, -edge/2.0, edge/2.0) 126 fdist += farthest*farthest 127 if fdist >= self.radius2: 128 return 2 # Box extends outside sphere 129 return 1 # Box entirely inside sphere 130 for i in range(3): 131 if not self.check_box[i]: continue 132 if right_edge[i] < self.bbox[i][0] or \ 133 left_edge[i] > self.bbox[i][1]: 134 return 0 # Box outside sphere bounding box 135 # http://www.gamedev.net/topic/335465-is-this-the-simplest-sphere-aabb-collision-test/ 136 cdist = 0 137 fdist = 0 138 for i in range(3): 139 # Early terminate 140 box_center = (right_edge[i] + left_edge[i])/2.0 141 relcenter = self.periodic_difference(box_center, self.center[i], i) 142 edge = right_edge[i] - left_edge[i] 143 closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0) 144 if relcenter >= 0: 145 farthest = relcenter + edge/2.0 146 else: 147 farthest = relcenter - edge/2.0 148 #farthest = relcenter + fclip(relcenter, -edge/2.0, edge/2.0) 149 cdist += closest*closest 150 fdist += farthest*farthest 151 if cdist > self.radius2: 152 return 0 # Box does not overlap sphere 153 if fdist < self.radius2: 154 return 1 # Sphere extends to entirely contain box 155 else: 156 return 2 # Sphere only partially overlaps box 157 158 def _hash_vals(self): 159 return (("radius", self.radius), 160 ("radius2", self.radius2), 161 ("center[0]", self.center[0]), 162 ("center[1]", self.center[1]), 163 ("center[2]", self.center[2])) 164 165 def _get_state_attnames(self): 166 return ("radius", "radius2", "center", "check_box") 167 168 def __setstate__(self, hashes): 169 super(SphereSelector, self).__setstate__(hashes) 170 self.set_bbox(self.center) 171 172 173sphere_selector = SphereSelector 174