Changeset 608
- Timestamp:
- Nov 22, 2004, 3:59:43 PM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/interpolate_sww.py
r600 r608 73 73 self.point_coordinates = self.read_xya(file_name) 74 74 self.interp.build_interpolation_matrix_A(self.point_coordinates) 75 #print "self.interp.A",self.interp.A76 75 77 76 #Change this if you want to generalise this application - -
inundation/ga/storm_surge/pyvolution/least_squares.py
r606 r608 239 239 #Find vertices near x 240 240 candidate_vertices = root.search(x[0], x[1]) 241 241 242 element_found, sigma0, sigma1, sigma2, k = \ 243 self.search_triangles_of_vertices(candidate_vertices, x) 244 245 is_more_elements = True 246 while not element_found and is_more_elements: 247 candidate_vertices = root.expand_search() 248 if candidate_vertices == []: 249 # All the triangles have been searched. 250 is_more_elements = False 251 else: 252 element_found, sigma0, sigma1, sigma2, k = \ 253 self.search_triangles_of_vertices(candidate_vertices, x) 254 255 242 256 257 #Update interpolation matrix A if necessary 258 if element_found is True: 259 #Assign values to matrix A 260 261 j0 = self.mesh.triangles[k,0] #Global vertex id 262 #self.A[i, j0] = sigma0 263 264 j1 = self.mesh.triangles[k,1] #Global vertex id 265 #self.A[i, j1] = sigma1 266 267 j2 = self.mesh.triangles[k,2] #Global vertex id 268 #self.A[i, j2] = sigma2 269 270 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} 271 js = [j0,j1,j2] 272 273 for j in js: 274 self.A[i,j] = sigmas[j] 275 for k in js: 276 self.AtA[j,k] += sigmas[j]*sigmas[k] 277 else: 278 pass 279 #Ok if there is no triangle for datapoint 280 #(as in brute force version) 281 #raise 'Could not find triangle for point', x 282 283 284 def search_triangles_of_vertices(self, candidate_vertices, x): 243 285 #Find triangle containing x: 244 element_found = False 245 286 element_found = False 287 288 # This will be returned if element_found = False 289 sigma2 = -10.0 290 sigma0 = -10.0 291 sigma1 = -10.0 292 246 293 #For all vertices in same cell as point x 247 294 for v in candidate_vertices: … … 255 302 xi2 = self.mesh.get_vertex_coordinate(k, 2) 256 303 304 #print "PDSG - k", k 305 #print "PDSG - xi0", xi0 306 #print "PDSG - xi1", xi1 307 #print "PDSG - xi2", xi2 308 #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))" \ 309 # % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1]) 310 257 311 #Get the three normals 258 312 n0 = self.mesh.get_normal(k, 0) … … 267 321 sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1) 268 322 323 #print "PDSG - sigma0", sigma0 324 #print "PDSG - sigma1", sigma1 325 #print "PDSG - sigma2", sigma2 326 269 327 #FIXME: Maybe move out to test or something 270 328 epsilon = 1.0e-6 … … 275 333 #Sigmas can get negative within 276 334 #machine precision on some machines (e.g nautilus) 277 #Hence the small eps 278 335 #Hence the small eps 279 336 eps = 1.0e-15 280 337 if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps: … … 285 342 #Don't look for any other triangle 286 343 break 344 return element_found, sigma0, sigma1, sigma2, k 287 345 288 289 290 291 #Update interpolation matrix A if necessary292 if element_found is True:293 #Assign values to matrix A294 295 j0 = self.mesh.triangles[k,0] #Global vertex id296 #self.A[i, j0] = sigma0297 298 j1 = self.mesh.triangles[k,1] #Global vertex id299 #self.A[i, j1] = sigma1300 301 j2 = self.mesh.triangles[k,2] #Global vertex id302 #self.A[i, j2] = sigma2303 304 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}305 js = [j0,j1,j2]306 307 for j in js:308 self.A[i,j] = sigmas[j]309 for k in js:310 self.AtA[j,k] += sigmas[j]*sigmas[k]311 else:312 pass313 #Ok if there is no triangle for datapoint314 #(as in brute force version)315 #raise 'Could not find triangle for point', x316 317 346 318 347 -
inundation/ga/storm_surge/pyvolution/quad.py
r484 r608 70 70 self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,self.name+'_se')) 71 71 self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,self.name+'_sw')) 72 72 73 73 74 75 74 def search(self, x, y): 75 #def search_new(self, x, y): 76 76 """Find all point indices sharing the same cell as point (x, y) 77 77 """ 78 78 branch = [] 79 79 points = [] 80 80 if self.children: 81 81 for child in self: 82 82 if child.contains(x,y): 83 points += child.search(x,y) 83 branch.append(self) 84 points, branch = child.search_branch(x,y, branch) 84 85 else: 85 86 # Leaf node: Get actual waypoints 86 87 points = self.retrieve() 87 88 89 self.branch = branch 88 90 return points 89 91 92 93 def search_branch(self, x, y, branch): 94 """Find all point indices sharing the same cell as point (x, y) 95 """ 96 points = [] 97 if self.children: 98 for child in self: 99 if child.contains(x,y): 100 branch.append(self) 101 points, branch = child.search_branch(x,y, branch) 102 103 else: 104 # Leaf node: Get actual waypoints 105 points = self.retrieve() 106 return points, branch 107 108 109 def expand_search(self): 110 """Find all point indices 'up' one cell from the last search 111 """ 112 points = [] 113 if self.branch == []: 114 points = [] 115 else: 116 larger_cell = self.branch.pop() 117 #print "larger_cell ", larger_cell.show() # Get_tree() 118 for child in larger_cell: 119 #This means it will also go down the branch known to be bad. 120 #So, if things are too slow, this could be speed up 121 points += child.retrieve() 122 return points 90 123 91 124 … … 367 400 #Hence, the root cell needs to be expanded slightly 368 401 ymax += (ymax-ymin)/10 369 xmax += (xmax-xmin)/10 402 xmax += (xmax-xmin)/10 370 403 371 404 -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r570 r608 40 40 41 41 42 def test_quad_tree(self): 43 p0 = [-10.0, -10.0] 44 p1 = [20.0, -10.0] 45 p2 = [-10.0, 20.0] 46 p3 = [10.0, 50.0] 47 p4 = [30.0, 30.0] 48 p5 = [50.0, 10.0] 49 p6 = [40.0, 60.0] 50 p7 = [60.0, 40.0] 51 p8 = [-66.0, 20.0] 52 p9 = [10.0, -66.0] 53 54 points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9] 55 triangles = [ [0, 1, 2], 56 [3, 2, 4], 57 [4, 2, 1], 58 [4, 1, 5], 59 [3, 4, 6], 60 [6, 4, 7], 61 [7, 4, 5], 62 [8, 0, 2], 63 [0, 9, 1]] 64 65 data = [ [4,4] ] 66 67 interp = Interpolation(points, triangles, data, alpha = 0.0) 68 #print "PDSG - interp.get_A()", interp.get_A() 69 answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0., 70 0., 0. , 0., 0., 0., 0.]] 71 assert allclose(interp.get_A(), answer) 72 42 73 43 74 def test_datapoints_at_vertices(self):
Note: See TracChangeset
for help on using the changeset viewer.