Changeset 6544
- Timestamp:
- Mar 18, 2009, 2:34:46 PM (16 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 5 edited
-
fit_interpolate/search_functions.py (modified) (7 diffs)
-
fit_interpolate/test_search_functions.py (modified) (10 diffs)
-
utilities/polygon.py (modified) (2 diffs)
-
utilities/polygon_ext.c (modified) (3 diffs)
-
utilities/test_polygon.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/search_functions.py
r6541 r6544 19 19 search_more_cells_time = initial_search_value 20 20 21 #FIXME test what happens if a 22 LAST_TRIANGLE = [[-10,[(num.array([max_float,max_float]),23 num.array([max_float,max_float]),24 num.array([max_float,max_float])),25 26 27 num.array([-1.1,-1.1]))]]]21 LAST_TRIANGLE = [[-10, 22 (num.array([[max_float,max_float], 23 [max_float,max_float], 24 [max_float,max_float]]), 25 (num.array([1,1], num.Int), #array default# 26 num.array([0,0], num.Int), #array default# 27 num.array([-1.1,-1.1])))]] 28 28 29 29 def search_tree_of_vertices(root, mesh, x): … … 48 48 global search_more_cells_time 49 49 50 #Find triangle containing x:51 element_found = False52 53 50 # This will be returned if element_found = False 51 element_found = False 54 52 sigma2 = -10.0 55 53 sigma0 = -10.0 … … 61 59 element_found, sigma0, sigma1, sigma2, k = \ 62 60 _search_triangles_of_vertices(mesh, last_triangle, x) 61 63 62 except: 63 print 'This should never happen:', last_triangle 64 64 element_found = False 65 65 66 #print "last_triangle", last_triangle67 66 if element_found is True: 68 #print "last_triangle", last_triangle69 67 return element_found, sigma0, sigma1, sigma2, k 70 68 71 # This was only slightly faster than just checking the 72 # last triangle and it significantly slowed down 73 # non-gridded fitting 74 # If the last element was a dud, search its neighbours 75 #print "last_triangle[0][0]", last_triangle[0][0] 76 #neighbours = mesh.get_triangle_neighbours(last_triangle[0][0]) 77 #print "neighbours", neighbours 78 #neighbours = [] 79 # for k in neighbours: 80 # if k >= 0: 81 # tri = mesh.get_vertex_coordinates(k, 82 # absolute=True) 83 # n0 = mesh.get_normal(k, 0) 84 # n1 = mesh.get_normal(k, 1) 85 # n2 = mesh.get_normal(k, 2) 86 # triangle =[[k,(tri, (n0, n1, n2))]] 87 # element_found, sigma0, sigma1, sigma2, k = \ 88 # _search_triangles_of_vertices(mesh, 89 # triangle, x) 90 # if element_found is True: 91 # return element_found, sigma0, sigma1, sigma2, k 92 93 #t0 = time.time() 69 94 70 # Get triangles in the cell that the point is in. 95 71 # Triangle is a list, first element triangle_id, … … 97 73 triangles = root.search(x[0], x[1]) 98 74 is_more_elements = True 99 100 75 element_found, sigma0, sigma1, sigma2, k = \ 101 76 _search_triangles_of_vertices(mesh, … … 133 108 global last_triangle 134 109 135 # these statments are needed if triangles is empty 136 #Find triangle containing x: 110 # These statments are needed if triangles is empty 137 111 element_found = False 138 139 # This will be returned if element_found = False140 112 sigma2 = -10.0 141 113 sigma0 = -10.0 142 114 sigma1 = -10.0 143 115 k = -10 144 145 # For all vertices in same cell as point x116 117 # For all vertices in same cell as point x 146 118 for k, tri_verts_norms in triangles: 147 119 tri = tri_verts_norms[0] … … 150 122 # tri is a list of verts (x, y), representing a tringle 151 123 # Find triangle that contains x (if any) and interpolate 152 if is_inside_triangle(x, tri, closed=True): 124 125 # Input check disabled to speed things up. 126 if is_inside_triangle(x, tri, 127 closed=True, 128 check_inputs=False): 129 153 130 n0, n1, n2 = tri_verts_norms[1] 154 element_found, sigma0, sigma1, sigma2 =\ 155 find_triangle_compute_interpolation(tri, n0, n1, n2, x) 156 else: 157 element_found = False 131 sigma0, sigma1, sigma2 =\ 132 compute_interpolation_values(tri, n0, n1, n2, x) 133 134 element_found = True 135 136 # Don't look for any other triangles in the triangle list 137 last_triangle = [[k, tri_verts_norms]] 138 break 158 139 159 140 160 if element_found is True:161 # Don't look for any other triangles in the triangle list162 last_triangle = [[k,tri_verts_norms]]163 break164 141 return element_found, sigma0, sigma1, sigma2, k 165 142 166 167 143 168 def find_triangle_compute_interpolation(triangle, n0, n1, n2, x): 169 """Compute linear interpolation of point x and triangle k in mesh. 170 It is assumed that x belongs to triangle k.max_float 144 def compute_interpolation_values(triangle, n0, n1, n2, x): 145 """Compute linear interpolation of point x and triangle. 146 147 n0, n1, n2 are normal to the tree edges. 171 148 """ 172 149 … … 174 151 xi0, xi1, xi2 = triangle 175 152 176 # this is where we can call some fast c code. 177 178 # Integrity check - machine precision is too hard 179 # so we use hardwired single precision 180 #epsilon = 1.0e-6 153 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0) 154 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1) 155 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2) 156 157 return sigma0, sigma1, sigma2 158 181 159 182 #if x[0] > max(xi0[0], xi1[0], xi2[0]) + epsilon:183 # return False,0,0,0184 #if x[0] < min(xi0[0], xi1[0], xi2[0]) - epsilon:185 # return False,0,0,0186 #if x[1] > max(xi0[1], xi1[1], xi2[1]) + epsilon:187 # return False,0,0,0188 #if x[1] < min(xi0[1], xi1[1], xi2[1]) - epsilon:189 # return False,0,0,0190 191 # machine precision on some machines (e.g. nautilus)192 epsilon = get_machine_precision() * 2193 194 # Compute interpolation - return as soon as possible195 # print "(xi0-xi1)", (xi0-xi1)196 # print "n0", n0197 # print "dot((xi0-xi1), n0)", dot((xi0-xi1), n0)198 199 sigma0 = num.dot((x-xi1), n0)/num.dot((xi0-xi1), n0)200 #if sigma0 < -epsilon:201 #print 'sigma0', sigma0202 #sigma0 = 0.0203 #return False,0,0,0204 sigma1 = num.dot((x-xi2), n1)/num.dot((xi1-xi2), n1)205 #if sigma1 < -epsilon:206 #print 'sigma1', sigma1207 #sigma1 = 0.0208 #return False,0,0,0209 sigma2 = num.dot((x-xi0), n2)/num.dot((xi2-xi0), n2)210 #if sigma2 < -epsilon:211 #print 'sigma2', sigma2212 #sigma2 = 0.0213 #return False,0,0,0214 215 # epsilon = 1.0e-6216 # we want to speed this up, so don't do assertions217 #delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero218 #msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\219 # %(delta, epsilon)220 #assert delta < epsilon, msg221 222 223 # Check that this triangle contains the data point224 # Sigmas are allowed to get negative within225 # machine precision on some machines (e.g. nautilus)226 #if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:227 # element_found = True228 #else:229 # element_found = False230 return True, sigma0, sigma1, sigma2231 232 160 def set_last_triangle(): 233 161 global last_triangle -
anuga_core/source/anuga/fit_interpolate/test_search_functions.py
r6541 r6544 5 5 from search_functions import search_tree_of_vertices, set_last_triangle 6 6 from search_functions import _search_triangles_of_vertices 7 from search_functions import find_triangle_compute_interpolation7 from search_functions import compute_interpolation_values 8 8 9 9 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh … … 11 11 from anuga.utilities.polygon import is_inside_polygon 12 12 from anuga.utilities.quad import build_quadtree, Cell 13 from anuga.utilities.numerical_tools import ensure_numeric 13 14 from anuga.utilities.polygon import is_inside_polygon, is_inside_triangle 14 15 … … 52 53 53 54 x = [0.2, 0.7] 54 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) 55 found, s0, s1, s2, k = search_tree_of_vertices(root, 56 mesh, 57 ensure_numeric(x)) 55 58 assert k == 1 # Triangle one 56 59 assert found is True 57 60 58 61 def test_bigger(self): 59 """test_larger mesh 62 """test_bigger 63 64 test larger mesh 60 65 """ 61 66 … … 72 77 [0.1,0.9], [0.4,0.6], [0.9,0.1], 73 78 [10, 3]]: 74 75 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) 79 80 81 found, s0, s1, s2, k = search_tree_of_vertices(root, 82 mesh, 83 ensure_numeric(x)) 76 84 77 85 if k >= 0: … … 95 103 mesh.check_integrity() 96 104 97 98 105 for m in range(8): 99 106 root = build_quadtree(mesh, max_points_per_cell = m) … … 105 112 [10, 3]]: 106 113 107 found, s0, s1, s2, k = search_tree_of_vertices(root, mesh, x) 114 found, s0, s1, s2, k = search_tree_of_vertices(root, 115 mesh, 116 x) 108 117 109 118 if k >= 0: 110 119 V = mesh.get_vertex_coordinates(k) # nodes for triangle k 120 121 assert is_inside_triangle(x, V, closed=True) 111 122 assert is_inside_polygon(x, V) 123 112 124 assert found is True 113 125 else: 114 126 assert found is False 115 127 116 128 129 if m == k == 0: return 117 130 118 131 def test_underlying_function(self): … … 127 140 128 141 # One point 129 x = [0.5, 0.5]142 x = ensure_numeric([0.5, 0.5]) 130 143 candidate_vertices = root.search(x[0], x[1]) 131 144 132 # print x, candidate_vertices145 # print x, candidate_vertices 133 146 found, sigma0, sigma1, sigma2, k = \ 134 147 _search_triangles_of_vertices(mesh, … … 156 169 _search_triangles_of_vertices(mesh, 157 170 candidate_vertices, 158 x)171 ensure_numeric(x)) 159 172 if k >= 0: 160 173 V = mesh.get_vertex_coordinates(k) # nodes for triangle k … … 211 224 212 225 213 def test_triangle_compute_interpolation(self): 214 """test_triangle_compute_interpolation 215 216 Test that triangle can be found if point is inside it 226 def test_compute_interpolation_values(self): 227 """test_compute_interpolation_values 228 229 Test that interpolation values are correc 230 231 This test used to check element_found as output from 232 find_triangle_compute_interpolation(triangle, n0, n1, n2, x) 233 and that failed before 18th March 2009. 234 235 Now this function no longer returns this flag, so the test 236 is merely checknig the sigmas. 237 217 238 """ 218 239 … … 232 253 closed=True, verbose=False) 233 254 assert is_inside_triangle(x, triangle, 234 closed=True, verbose=False) 235 236 element_found,sigma0, sigma1, sigma2 = \237 find_triangle_compute_interpolation(triangle, n0, n1, n2, x)255 closed=True, verbose=False) 256 257 sigma0, sigma1, sigma2 = \ 258 compute_interpolation_values(triangle, n0, n1, n2, x) 238 259 239 260 msg = 'Point which is clearly inside triangle was not found' 240 assert element_found is True, msg 241 242 #print sigma0, sigma1, sigma2 243 #print sigma0 + sigma1 + sigma2 244 245 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6 261 assert abs(sigma0 + sigma1 + sigma2 - 1) < 1.0e-6, msg 246 262 247 263 #------------------------------------------------------------- -
anuga_core/source/anuga/utilities/polygon.py
r6541 r6544 232 232 """ 233 233 234 triangle = ensure_numeric(triangle) 235 point = ensure_numeric(point, num.Float) 236 234 237 if check_inputs is True: 235 triangle = ensure_numeric(triangle)236 237 point = ensure_numeric(point, num.Float) # Convert point to array of points238 238 msg = 'is_inside_triangle must be invoked with one point only' 239 239 assert num.allclose(point.shape, [2]), msg 240 240 241 241 242 # Use C-implementation 243 return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol)) 244 245 246 247 # FIXME (Ole): The rest of this function has been made 248 # obsolete by the C extension. 249 242 250 # Quickly reject points that are clearly outside 243 251 if point[0] < min(triangle[:,0]): return False … … 245 253 if point[1] < min(triangle[:,1]): return False 246 254 if point[1] > max(triangle[:,1]): return False 247 248 # Use C-implementation249 return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))250 251 252 253 # FIXME (Ole): The rest of this function has been made254 # obsolete by the C extension.255 255 256 256 # Start search -
anuga_core/source/anuga/utilities/polygon_ext.c
r6535 r6544 254 254 double a00, a10, a01, a11, b0, b1; 255 255 double denom, alpha, beta; 256 257 double x, y; // Point coordinates 256 258 int i, j, res; 257 259 260 x = point[0]; 261 y = point[1]; 262 263 // Quickly reject points that are clearly outside 264 if ((x < triangle[0]) && 265 (x < triangle[2]) && 266 (x < triangle[4])) return 0; 267 268 if ((x > triangle[0]) && 269 (x > triangle[2]) && 270 (x > triangle[4])) return 0; 271 272 if ((y < triangle[1]) && 273 (y < triangle[3]) && 274 (y < triangle[5])) return 0; 275 276 if ((y > triangle[1]) && 277 (y > triangle[3]) && 278 (y > triangle[5])) return 0; 279 280 258 281 // v0 = C-A 259 282 v0x = triangle[4]-triangle[0]; … … 274 297 if (fabs(denom) > 0.0) { 275 298 // v = point-A 276 vx = point[0]- triangle[0];277 vy = point[1]- triangle[1];299 vx = x - triangle[0]; 300 vy = y - triangle[1]; 278 301 279 302 b0 = v0x*vx + v0y*vy; // innerproduct(v0, v) … … 291 314 for (i=0; i<3; i++) { 292 315 j = (i+1) % 3; // Circular index into triangle array 293 res = __point_on_line( point[0], point[1],316 res = __point_on_line(x, y, 294 317 triangle[2*i], triangle[2*i+1], 295 318 triangle[2*j], triangle[2*j+1], -
anuga_core/source/anuga/utilities/test_polygon.py
r6535 r6544 1896 1896 assert res is False 1897 1897 1898 1898 1899 res = is_inside_triangle([10, 3], 1900 [[0.1, 0.], 1901 [0.1, 0.08333333], 1902 [0., 0.]]) 1903 assert res is False 1904 1899 1905 1900 1906 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
