Changeset 6410
- Timestamp:
- Feb 25, 2009, 9:37:22 AM (15 years ago)
- Location:
- branches/numpy/anuga
- Files:
-
- 47 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py
r6304 r6410 377 377 #indices = range(M) 378 378 379 return num.take(self.triangles, indices )379 return num.take(self.triangles, indices, axis=0) 380 380 381 381 -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py
r6304 r6410 479 479 indices, verbose=verbose, 480 480 use_cache=use_cache) 481 # if (isinstance(numeric, float) or isinstance(numeric, int)482 # or isinstance(numeric, long)):483 481 else: 484 482 try: 485 483 numeric = float(numeric) 486 484 except ValueError: 487 msg = ("Illegal type for variable 'numeric': " 488 "%s, shape=%s\n%s" 489 % (type(numeric), numeric.shape, str(numeric))) 490 msg += ('\ntype(numeric)==FloatType is %s' 491 % str(type(numeric)==FloatType)) 492 msg += ('\nisinstance(numeric, float)=%s' 493 % str(isinstance(numeric, float))) 494 msg += ('\nisinstance(numeric, num.ndarray)=%s' 495 % str(isinstance(numeric, num.ndarray))) 485 msg = ("Illegal type for variable 'numeric': %s" 486 % type(numeric)) 496 487 raise Exception, msg 497 488 self.set_values_from_constant(numeric, location, 498 489 indices, verbose) 499 elif quantity is not None:500 if type(numeric) in [FloatType, IntType, LongType]:501 self.set_values_from_constant(numeric, location,502 indices, verbose)503 else:504 msg = ("Illegal type for variable 'numeric': %s, shape=%s\n%s"505 % (type(numeric), numeric.shape, str(numeric)))506 msg += ('\ntype(numeric)==FloatType is %s'507 % str(type(numeric)==FloatType))508 msg += ('\nisinstance(numeric, float)=%s'509 % str(isinstance(numeric, float)))510 msg += ('\nisinstance(numeric, num.ndarray)=%s'511 % str(isinstance(numeric, num.ndarray)))512 raise Exception, msg513 490 elif quantity is not None: 514 491 self.set_values_from_quantity(quantity, location, indices, verbose) … … 745 722 indices = range(len(self)) 746 723 747 V = num.take(self.domain.get_centroid_coordinates(), indices )724 V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0) 748 725 x = V[:,0]; y = V[:,1] 749 726 if use_cache is True: … … 1213 1190 if (indices == None): 1214 1191 indices = range(len(self)) 1215 return num.take(self.centroid_values, indices )1192 return num.take(self.centroid_values, indices, axis=0) 1216 1193 elif location == 'edges': 1217 1194 if (indices == None): 1218 1195 indices = range(len(self)) 1219 return num.take(self.edge_values, indices )1196 return num.take(self.edge_values, indices, axis=0) 1220 1197 elif location == 'unique vertices': 1221 1198 if (indices == None): … … 1243 1220 if (indices is None): 1244 1221 indices = range(len(self)) 1245 return num.take(self.vertex_values, indices )1222 return num.take(self.vertex_values, indices, axis=0) 1246 1223 1247 1224 ## -
branches/numpy/anuga/abstract_2d_finite_volumes/quantity_ext.c
r6304 r6410 1030 1030 } 1031 1031 1032 // check that numpy array objects arrays are C contiguous memory 1033 CHECK_C_CONTIG(vertex_value_indices); 1034 CHECK_C_CONTIG(number_of_triangles_per_node); 1035 CHECK_C_CONTIG(vertex_values); 1036 CHECK_C_CONTIG(A); 1037 1032 1038 N = vertex_value_indices -> dimensions[0]; 1033 1039 // printf("Got parameters, N=%d\n", N); -
branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py
r6304 r6410 764 764 765 765 #print domain.quantities['friction'].get_values() 766 msg = ("domain.quantities['friction'].get_values()=\n%s\n" 767 'should equal\n' 768 '[[ 0.09, 0.09, 0.09],\n' 769 ' [ 0.09, 0.09, 0.09],\n' 770 ' [ 0.07, 0.07, 0.07],\n' 771 ' [ 0.07, 0.07, 0.07],\n' 772 ' [ 1.0, 1.0, 1.0],\n' 773 ' [ 1.0, 1.0, 1.0]]' 774 % str(domain.quantities['friction'].get_values())) 766 775 assert num.allclose(domain.quantities['friction'].get_values(), 767 776 [[ 0.09, 0.09, 0.09], … … 770 779 [ 0.07, 0.07, 0.07], 771 780 [ 1.0, 1.0, 1.0], 772 [ 1.0, 1.0, 1.0]]) 781 [ 1.0, 1.0, 1.0]]), msg 773 782 774 783 domain.set_region([set_bottom_friction, set_top_friction]) … … 793 802 794 803 #------------------------------------------------------------- 804 795 805 if __name__ == "__main__": 796 806 suite = unittest.makeSuite(Test_Domain,'test') -
branches/numpy/anuga/abstract_2d_finite_volumes/test_ermapper.py
r6304 r6410 183 183 184 184 #------------------------------------------------------------- 185 185 186 if __name__ == "__main__": 186 187 suite = unittest.makeSuite(Test_ERMapper,'test') 187 188 runner = unittest.TextTestRunner() 188 189 runner.run(suite) 189 190 191 192 -
branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py
r6304 r6410 63 63 domain = General_mesh(nodes, triangles, 64 64 geo_reference = geo) 65 65 66 verts = domain.get_vertex_coordinates(triangle_id=0) 66 self.assert_(num.allclose(num.array([b,a,c]), verts)) 67 msg = ("num.array([b,a,c])=\n%s\nshould be the same as 'verts'=\n%s" 68 % (str(num.array([b,a,c])), str(verts))) 69 #self.assert_(num.allclose(num.array([b,a,c]), verts)) 70 assert num.allclose(num.array([b,a,c]), verts), msg 71 67 72 verts = domain.get_vertex_coordinates(triangle_id=0) 68 73 self.assert_(num.allclose(num.array([b,a,c]), verts)) … … 117 122 domain = General_mesh(nodes, triangles) 118 123 119 value = [7] 120 assert num.allclose(domain.get_triangles(), triangles) 124 # value = [7] 125 msg = ("domain.get_triangles()=\n%s\nshould be the same as " 126 "'triangles'=\n%s" 127 % (str(domain.get_triangles()), str(triangles))) 128 assert num.allclose(domain.get_triangles(), triangles), msg 129 msg = ("domain.get_triangles([0,4])=\n%s\nshould be the same as " 130 "'[triangles[0], triangles[4]]' which is\n%s" 131 % (str(domain.get_triangles([0,4])), 132 str([triangles[0], triangles[4]]))) 121 133 assert num.allclose(domain.get_triangles([0,4]), 122 [triangles[0], triangles[4]])134 [triangles[0], triangles[4]]), msg 123 135 124 136 … … 278 290 node = domain.get_node(2) 279 291 #self.assertEqual(c, node) 280 self.failUnless(num.alltrue(c == node)) 292 msg = ('\nc=%s\nmode=%s' % (str(c), str(node))) 293 self.failUnless(num.alltrue(c == node), msg) 281 294 282 295 node = domain.get_node(2, absolute=True) … … 328 341 329 342 #------------------------------------------------------------- 343 330 344 if __name__ == "__main__": 331 345 suite = unittest.makeSuite(Test_General_Mesh,'test') 332 #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')333 346 runner = unittest.TextTestRunner() 334 347 runner.run(suite) -
branches/numpy/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
r6304 r6410 270 270 271 271 272 273 274 275 272 #------------------------------------------------------------- 273 276 274 if __name__ == "__main__": 277 275 suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test') 278 276 runner = unittest.TextTestRunner() 279 277 runner.run(suite) 280 281 282 -
branches/numpy/anuga/abstract_2d_finite_volumes/test_ghost.py
r6304 r6410 45 45 46 46 47 #------------------------------------------------------------- 47 48 48 49 #-------------------------------------------------------------50 49 if __name__ == "__main__": 51 50 suite = unittest.makeSuite(Test_Domain,'test') -
branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
r6304 r6410 1839 1839 1840 1840 1841 1842 1843 1841 #------------------------------------------------------------- 1842 1844 1843 if __name__ == "__main__": 1845 #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')1846 #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')1847 #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')1848 1844 suite = unittest.makeSuite(Test_Mesh,'test') 1849 1845 runner = unittest.TextTestRunner() 1850 1846 runner.run(suite) 1851 1852 1853 1854 -
branches/numpy/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py
r6304 r6410 255 255 256 256 #------------------------------------------------------------- 257 257 258 if __name__ == "__main__": 258 259 suite = unittest.makeSuite(Test_pmesh2domain, 'test') 259 260 runner = unittest.TextTestRunner() 260 261 runner.run(suite) 261 262 263 264 -
branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py
r6304 r6410 504 504 quantity.set_values(0.0) 505 505 quantity.set_values(3.14, polygon=polygon) 506 msg = ('quantity.vertex_values=\n%s\nshould be close to\n' 507 '[[0,0,0],\n' 508 ' [3.14,3.14,3.14],\n' 509 ' [3.14,3.14,3.14],\n' 510 ' [0,0,0]]' % str(quantity.vertex_values)) 506 511 assert num.allclose(quantity.vertex_values, 507 512 [[0,0,0], 508 513 [3.14,3.14,3.14], 509 514 [3.14,3.14,3.14], 510 [0,0,0]]) 515 [0,0,0]]), msg 511 516 512 517 … … 2508 2513 2509 2514 #------------------------------------------------------------- 2515 2510 2516 if __name__ == "__main__": 2511 2517 suite = unittest.makeSuite(Test_Quantity, 'test') 2512 #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')2513 2518 runner = unittest.TextTestRunner() 2514 2519 runner.run(suite) -
branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py
r6304 r6410 189 189 190 190 def test_unique_vertices_average_loc_vert(self): 191 """ 192 get values based on triangle lists. 193 """194 from mesh_factory import rectangular195 from shallow_water import Domain 196 197 #Create basic mesh198 points, vertices, boundary = rectangular(1, 3) 199 #Create shallow water domain 200 domain = Domain(points, vertices, boundary) 201 domain.build_tagged_elements_dictionary({'bottom': [0,1],202 'top': [4,5],203 'not_bottom': [2,3,4,5]})191 """Get values based on triangle lists.""" 192 193 from mesh_factory import rectangular 194 from shallow_water import Domain 195 196 #Create basic mesh 197 points, vertices, boundary = rectangular(1, 3) 198 199 #Create shallow water domain 200 domain = Domain(points, vertices, boundary) 201 domain.build_tagged_elements_dictionary({'bottom': [0, 1], 202 'top': [4, 5], 203 'not_bottom': [2, 3, 4, 5]}) 204 204 205 205 #Set friction 206 206 domain.set_quantity('friction', add_x_y) 207 av_bottom = 2.0 /3.0207 av_bottom = 2.0 / 3.0 208 208 add = 60.0 209 209 calc_frict = av_bottom + add 210 210 domain.set_region(Add_value_to_region('bottom', 'friction', add, 211 initial_quantity='friction', 212 #location='unique vertices', 213 location='vertices', 214 average=True 215 )) 216 217 #print domain.quantities['friction'].get_values() 211 initial_quantity='friction', 212 location='vertices', 213 average=True)) 214 218 215 frict_points = domain.quantities['friction'].get_values() 216 219 217 expected = [calc_frict, calc_frict, calc_frict] 220 msg = ('frict_points[0]=%s\nexpected=%s' % (str(frict_points[0]),221 218 msg = ('frict_points[0]=%s\nexpected=%s' 219 % (str(frict_points[0]), str(expected))) 222 220 assert num.allclose(frict_points[0], expected), msg 223 msg = ('frict_points[1]=%s\nexpected=%s' % (str(frict_points[1]), 224 str(expected))) 221 222 msg = ('frict_points[1]=%s\nexpected=%s' 223 % (str(frict_points[1]), str(expected))) 225 224 assert num.allclose(frict_points[1], expected), msg 226 225 … … 264 263 265 264 #------------------------------------------------------------- 265 266 266 if __name__ == "__main__": 267 suite = unittest.makeSuite(Test_Region, 'test') 267 #suite = unittest.makeSuite(Test_Region, 'test') 268 suite = unittest.makeSuite(Test_Region, 'test_unique_vertices_average_loc_vert') 268 269 runner = unittest.TextTestRunner() 269 270 runner.run(suite) -
branches/numpy/anuga/abstract_2d_finite_volumes/test_util.py
r6304 r6410 183 183 184 184 last_time_index = len(time)-1 #Last last_time_index 185 d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1)) 186 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1)) 187 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1)) 185 d_stage = num.reshape(num.take(stage[last_time_index, :], 186 [0,5,10,15], 187 axis=0), 188 (4,1)) 189 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], 190 [0,5,10,15], 191 axis=0), 192 (4,1)) 193 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], 194 [0,5,10,15], 195 axis=0), 196 (4,1)) 188 197 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 189 198 … … 195 204 196 205 #And the midpoints are found now 197 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15] )198 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15] )206 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0) 207 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0) 199 208 200 209 diag = num.concatenate( (Dx, Dy), axis=1) … … 219 228 220 229 timestep = 0 #First timestep 221 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))222 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))223 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))230 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 231 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 232 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 224 233 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 225 234 … … 241 250 242 251 timestep = 33 243 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))244 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))245 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))252 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 253 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 254 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 246 255 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 247 256 … … 262 271 263 272 timestep = 15 264 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))265 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))266 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))273 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 274 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 275 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 267 276 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 268 277 … … 275 284 # 276 285 timestep = 16 277 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15] ), (4,1))278 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15] ), (4,1))279 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15] ), (4,1))286 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1)) 287 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 288 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1)) 280 289 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 281 290 … … 386 395 387 396 last_time_index = len(time)-1 #Last last_time_index 388 d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1)) 389 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1)) 390 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1)) 397 d_stage = num.reshape(num.take(stage[last_time_index, :], 398 [0,5,10,15], 399 axis=0), 400 (4,1)) 401 d_uh = num.reshape(num.take(xmomentum[last_time_index, :], 402 [0,5,10,15], 403 axis=0), 404 (4,1)) 405 d_vh = num.reshape(num.take(ymomentum[last_time_index, :], 406 [0,5,10,15], 407 axis=0), 408 (4,1)) 391 409 D = num.concatenate((d_stage, d_uh, d_vh), axis=1) 392 410 … … 398 416 399 417 #And the midpoints are found now 400 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15] )401 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15] )418 Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0) 419 Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0) 402 420 403 421 diag = num.concatenate( (Dx, Dy), axis=1) … … 415 433 416 434 t = time[last_time_index] 417 q = f(t, point_id=0); assert num.allclose(r0, q) 418 q = f(t, point_id=1); assert num.allclose(r1, q) 419 q = f(t, point_id=2); assert num.allclose(r2, q) 435 436 q = f(t, point_id=0) 437 msg = '\nr0=%s\nq=%s' % (str(r0), str(q)) 438 assert num.allclose(r0, q), msg 439 440 q = f(t, point_id=1) 441 msg = '\nr1=%s\nq=%s' % (str(r1), str(q)) 442 assert num.allclose(r1, q), msg 443 444 q = f(t, point_id=2) 445 msg = '\nr2=%s\nq=%s' % (str(r2), str(q)) 446 assert num.allclose(r2, q), msg 420 447 421 448 … … 424 451 425 452 timestep = 0 #First timestep 426 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 427 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 428 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 453 d_stage = num.reshape(num.take(stage[timestep, :], 454 [0,5,10,15], 455 axis=0), 456 (4,1)) 457 d_uh = num.reshape(num.take(xmomentum[timestep, :], 458 [0,5,10,15], 459 axis=0), 460 (4,1)) 461 d_vh = num.reshape(num.take(ymomentum[timestep, :], 462 [0,5,10,15], 463 axis=0), 464 (4,1)) 429 465 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 430 466 … … 446 482 447 483 timestep = 33 448 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 449 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 450 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 484 d_stage = num.reshape(num.take(stage[timestep, :], 485 [0,5,10,15], 486 axis=0), 487 (4,1)) 488 d_uh = num.reshape(num.take(xmomentum[timestep, :], 489 [0,5,10,15], 490 axis=0), 491 (4,1)) 492 d_vh = num.reshape(num.take(ymomentum[timestep, :], 493 [0,5,10,15], 494 axis=0), 495 (4,1)) 451 496 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 452 497 … … 467 512 468 513 timestep = 15 469 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 470 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 471 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 514 d_stage = num.reshape(num.take(stage[timestep, :], 515 [0,5,10,15], 516 axis=0), 517 (4,1)) 518 d_uh = num.reshape(num.take(xmomentum[timestep, :], 519 [0,5,10,15], 520 axis=0), 521 (4,1)) 522 d_vh = num.reshape(num.take(ymomentum[timestep, :], 523 [0,5,10,15], 524 axis=0), 525 (4,1)) 472 526 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 473 527 … … 480 534 # 481 535 timestep = 16 482 d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1)) 483 d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1)) 484 d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1)) 536 d_stage = num.reshape(num.take(stage[timestep, :], 537 [0,5,10,15], 538 axis=0), 539 (4,1)) 540 d_uh = num.reshape(num.take(xmomentum[timestep, :], 541 [0,5,10,15], 542 axis=0), 543 (4,1)) 544 d_vh = num.reshape(num.take(ymomentum[timestep, :], 545 [0,5,10,15], 546 axis=0), 547 (4,1)) 485 548 D = num.concatenate( (d_stage, d_uh, d_vh), axis=1) 486 549 … … 1913 1976 if 314 < angle < 316: v=1 1914 1977 assert v==1 1915 1916 1917 1918 1919 1978 1920 1979 1921 1980 #------------------------------------------------------------- 1981 1922 1982 if __name__ == "__main__": 1923 1983 suite = unittest.makeSuite(Test_Util,'test') 1924 # suite = unittest.makeSuite(Test_Util,'test_sww2csv_gauges')1925 1984 # runner = unittest.TextTestRunner(verbosity=2) 1926 1985 runner = unittest.TextTestRunner(verbosity=1) 1927 1986 runner.run(suite) 1928 1929 1930 1931 -
branches/numpy/anuga/abstract_2d_finite_volumes/util.py
r6304 r6410 476 476 if boundary_polygon is not None: 477 477 #removes sts points that do not lie on boundary 478 quantities[name] = num.take(quantities[name], gauge_id, 1) 478 #quantities[name] = num.take(quantities[name], gauge_id, 1) #was# 479 quantities[name] = num.take(quantities[name], gauge_id, axis=1) 479 480 480 481 # Close sww, tms or sts netcdf file … … 1798 1799 # Remove the loners from verts 1799 1800 # Could've used X=compress(less(loners,N),loners) 1800 # verts=num.take(verts,X ) to Remove the loners from verts1801 # verts=num.take(verts,X,axis=0) to Remove the loners from verts 1801 1802 # but I think it would use more memory 1802 1803 new_i = lone_start # point at first loner - 'shuffle down' target -
branches/numpy/anuga/alpha_shape/test_alpha_shape.py
r6304 r6410 394 394 (46, 45)] 395 395 assert num.allclose(answer, result) 396 396 397 #------------------------------------------------------------- 398 397 399 if __name__ == "__main__": 398 #print "starting tests \n"399 400 suite = unittest.makeSuite(TestCase,'test') 400 401 runner = unittest.TextTestRunner(verbosity=1) 401 402 runner.run(suite) 402 #print "finished tests \n"403 404 405 406 407 408 409 -
branches/numpy/anuga/caching/test_caching.py
r6304 r6410 891 891 #------------------------------------------------------------- 892 892 if __name__ == "__main__": 893 #suite = unittest.makeSuite(Test_Caching, 'test_caching_of_callable_objects')894 893 suite = unittest.makeSuite(Test_Caching, 'test') 895 894 runner = unittest.TextTestRunner() -
branches/numpy/anuga/coordinate_transforms/test_geo_reference.py
r6360 r6410 505 505 506 506 #------------------------------------------------------------- 507 507 508 if __name__ == "__main__": 508 509 509 suite = unittest.makeSuite(geo_referenceTestCase,'test') 510 #suite = unittest.makeSuite(geo_referenceTestCase,'test_get_absolute')511 510 runner = unittest.TextTestRunner() #verbosity=2) 512 511 runner.run(suite) -
branches/numpy/anuga/coordinate_transforms/test_lat_long_UTM_conversion.py
r6304 r6410 121 121 assert num.allclose(lat, lat_calced) 122 122 assert num.allclose(lon, long_calced) 123 123 124 #------------------------------------------------------------- 125 124 126 if __name__ == "__main__": 125 126 #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')127 127 mysuite = unittest.makeSuite(TestCase,'test') 128 128 runner = unittest.TextTestRunner() -
branches/numpy/anuga/coordinate_transforms/test_redfearn.py
r6304 r6410 411 411 412 412 #------------------------------------------------------------- 413 413 414 if __name__ == "__main__": 414 415 #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')416 #mysuite = unittest.makeSuite(TestCase,'test_UTM_6_nonstandard_projection')417 415 mysuite = unittest.makeSuite(TestCase,'test') 418 416 runner = unittest.TextTestRunner() -
branches/numpy/anuga/fit_interpolate/interpolate.py
r6304 r6410 68 68 max_vertices_per_cell=None, 69 69 start_blocking_len=500000, 70 use_cache=False, 71 verbose=False): 70 use_cache=False, 71 verbose=False): 72 72 """Interpolate vertex_values to interpolation points. 73 73 74 74 Inputs (mandatory): 75 75 76 76 77 77 vertex_coordinates: List of coordinate pairs [xi, eta] of 78 points constituting a mesh 78 points constituting a mesh 79 79 (or an m x 2 numeric array or 80 80 a geospatial object) … … 83 83 84 84 triangles: List of 3-tuples (or a numeric array) of 85 integers representing indices of all vertices 85 integers representing indices of all vertices 86 86 in the mesh. 87 87 88 88 vertex_values: Vector or array of data at the mesh vertices. 89 89 If array, interpolation will be done for each column as 90 90 per underlying matrix-matrix multiplication 91 91 92 92 interpolation_points: Interpolate mesh data to these positions. 93 93 List of coordinate pairs [x, y] of 94 94 data points or an nx2 numeric array or a 95 95 Geospatial_data object 96 96 97 97 Inputs (optional) 98 98 99 99 mesh_origin: A geo_reference object or 3-tuples consisting of 100 100 UTM zone, easting and northing. … … 108 108 object and a mesh origin, since geospatial has its 109 109 own mesh origin. 110 110 111 111 start_blocking_len: If the # of points is more or greater than this, 112 start blocking 113 112 start blocking 113 114 114 use_cache: True or False 115 115 116 116 117 117 Output: 118 118 119 119 Interpolated values at specified point_coordinates 120 120 121 Note: This function is a simple shortcut for case where 121 Note: This function is a simple shortcut for case where 122 122 interpolation matrix is unnecessary 123 Note: This function does not take blocking into account, 123 Note: This function does not take blocking into account, 124 124 but allows caching. 125 125 126 126 """ 127 127 128 128 # FIXME(Ole): Probably obsolete since I is precomputed and 129 129 # interpolate_block caches 130 130 131 131 from anuga.caching import cache 132 132 133 133 # Create interpolation object with matrix 134 args = (ensure_numeric(vertex_coordinates, num.float), 134 args = (ensure_numeric(vertex_coordinates, num.float), 135 135 ensure_numeric(triangles)) 136 136 kwargs = {'mesh_origin': mesh_origin, 137 137 'max_vertices_per_cell': max_vertices_per_cell, 138 138 'verbose': verbose} 139 139 140 140 if use_cache is True: 141 141 if sys.platform != 'win32': … … 149 149 else: 150 150 I = apply(Interpolate, args, kwargs) 151 151 152 152 # Call interpolate method with interpolation points 153 153 result = I.interpolate_block(vertex_values, interpolation_points, 154 154 use_cache=use_cache, 155 155 verbose=verbose) 156 156 157 157 return result 158 158 159 159 160 160 ## 161 # @brief 161 # @brief 162 162 class Interpolate (FitInterpolate): 163 163 … … 181 181 Inputs: 182 182 vertex_coordinates: List of coordinate pairs [xi, eta] of 183 183 points constituting a mesh (or an m x 2 numeric array or 184 184 a geospatial object) 185 185 Points may appear multiple times … … 202 202 203 203 # FIXME (Ole): Need an input check 204 204 205 205 FitInterpolate.__init__(self, 206 206 vertex_coordinates=vertex_coordinates, … … 209 209 verbose=verbose, 210 210 max_vertices_per_cell=max_vertices_per_cell) 211 211 212 212 # Initialise variables 213 213 self._A_can_be_reused = False # FIXME (Ole): Probably obsolete … … 215 215 self.interpolation_matrices = {} # Store precomputed matrices 216 216 217 217 218 218 219 219 ## … … 243 243 point_coordinates: Interpolate mesh data to these positions. 244 244 List of coordinate pairs [x, y] of 245 246 247 245 data points or an nx2 numeric array or a Geospatial_data object 246 247 If point_coordinates is absent, the points inputted last time 248 248 this method was called are used, if possible. 249 249 250 250 start_blocking_len: If the # of points is more or greater than this, 251 start blocking 252 253 254 251 start blocking 252 253 Output: 254 Interpolated values at inputted points (z). 255 255 """ 256 256 … … 262 262 # This has now been addressed through an attempt in interpolate_block 263 263 264 if verbose: print 'Build intepolation object' 264 if verbose: print 'Build intepolation object' 265 265 if isinstance(point_coordinates, Geospatial_data): 266 266 point_coordinates = point_coordinates.get_data_points(absolute=True) … … 280 280 msg = 'ERROR (interpolate.py): No point_coordinates inputted' 281 281 raise Exception(msg) 282 282 283 283 if point_coordinates is not None: 284 284 self._point_coordinates = point_coordinates … … 293 293 start = 0 294 294 # creating a dummy array to concatenate to. 295 295 296 296 f = ensure_numeric(f, num.float) 297 297 if len(f.shape) > 1: … … 299 299 else: 300 300 z = num.zeros((0,), num.int) #array default# 301 301 302 302 for end in range(start_blocking_len, 303 303 len(point_coordinates), … … 313 313 z = num.concatenate((z, t)) 314 314 return z 315 315 316 316 317 317 ## … … 322 322 # @param verbose True if this function is verbose. 323 323 # @return ?? 324 def interpolate_block(self, f, point_coordinates, 324 def interpolate_block(self, f, point_coordinates, 325 325 use_cache=False, verbose=False): 326 326 """ … … 329 329 330 330 Return the point data, z. 331 331 332 332 See interpolate for doc info. 333 333 """ 334 334 335 335 # FIXME (Ole): I reckon we should change the interface so that 336 # the user can specify the interpolation matrix instead of the 336 # the user can specify the interpolation matrix instead of the 337 337 # interpolation points to save time. 338 338 … … 342 342 # Convert lists to numeric arrays if necessary 343 343 point_coordinates = ensure_numeric(point_coordinates, num.float) 344 f = ensure_numeric(f, num.float) 344 f = ensure_numeric(f, num.float) 345 345 346 346 from anuga.caching import myhash 347 import sys 348 347 import sys 348 349 349 if use_cache is True: 350 350 if sys.platform != 'win32': 351 351 # FIXME (Ole): (Why doesn't this work on windoze?) 352 352 # Still absolutely fails on Win 24 Oct 2008 353 353 354 354 X = cache(self._build_interpolation_matrix_A, 355 355 args=(point_coordinates,), 356 kwargs={'verbose': verbose}, 357 verbose=verbose) 356 kwargs={'verbose': verbose}, 357 verbose=verbose) 358 358 else: 359 359 # FIXME … … 361 361 # This will work on Linux as well if we want to use it there. 362 362 key = myhash(point_coordinates) 363 364 reuse_A = False 365 363 364 reuse_A = False 365 366 366 if self.interpolation_matrices.has_key(key): 367 X, stored_points = self.interpolation_matrices[key] 367 X, stored_points = self.interpolation_matrices[key] 368 368 if num.alltrue(stored_points == point_coordinates): 369 reuse_A = True 370 369 reuse_A = True # Reuse interpolation matrix 370 371 371 if reuse_A is False: 372 372 X = self._build_interpolation_matrix_A(point_coordinates, … … 375 375 else: 376 376 X = self._build_interpolation_matrix_A(point_coordinates, 377 verbose=verbose) 378 379 # Unpack result 377 verbose=verbose) 378 379 # Unpack result 380 380 self._A, self.inside_poly_indices, self.outside_poly_indices = X 381 381 … … 389 389 msg += ' I got %d points and %d matrix rows.' \ 390 390 % (point_coordinates.shape[0], self._A.shape[0]) 391 assert point_coordinates.shape[0] == self._A.shape[0], msg 391 assert point_coordinates.shape[0] == self._A.shape[0], msg 392 392 393 393 msg = 'The number of columns in matrix A must be the same as the ' 394 394 msg += 'number of mesh vertices.' 395 395 msg += ' I got %d vertices and %d matrix columns.' \ 396 % (f.shape[0], self._A.shape[1]) 396 % (f.shape[0], self._A.shape[1]) 397 397 assert self._A.shape[1] == f.shape[0], msg 398 398 … … 409 409 """ 410 410 Return the point data, z. 411 411 412 412 Precondition: The _A matrix has been created 413 413 """ … … 416 416 417 417 # Taking into account points outside the mesh. 418 for i in self.outside_poly_indices: 418 for i in self.outside_poly_indices: 419 419 z[i] = NAN 420 420 return z … … 456 456 self.mesh.get_boundary_polygon(), 457 457 closed=True, verbose=verbose) 458 458 459 459 # Build n x m interpolation matrix 460 460 if verbose and len(outside_poly_indices) > 0: 461 461 print '\n WARNING: Points outside mesh boundary. \n' 462 462 463 463 # Since you can block, throw a warning, not an error. 464 464 if verbose and 0 == len(inside_poly_indices): 465 465 print '\n WARNING: No points within the mesh! \n' 466 466 467 467 m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) 468 468 n = point_coordinates.shape[0] # Nbr of data points … … 474 474 475 475 n = len(inside_poly_indices) 476 476 477 477 # Compute matrix elements for points inside the mesh 478 478 if verbose: print 'Building interpolation matrix from %d points' %n … … 485 485 element_found, sigma0, sigma1, sigma2, k = \ 486 486 search_tree_of_vertices(self.root, self.mesh, x) 487 488 487 488 # Update interpolation matrix A if necessary 489 489 if element_found is True: 490 490 # Assign values to matrix A … … 499 499 A[i, j] = sigmas[j] 500 500 else: 501 msg = 'Could not find triangle for point', x 501 msg = 'Could not find triangle for point', x 502 502 raise Exception(msg) 503 503 return A, inside_poly_indices, outside_poly_indices 504 504 505 505 506 507 508 509 506 507 508 509 510 510 ## 511 511 # @brief ?? … … 527 527 List of coordinate pairs [x, y] of 528 528 data points or an nx2 numeric array or a Geospatial_data object 529 529 530 530 No test for this yet. 531 531 Note, this has no time the input data has no time dimension. Which is 532 532 different from most of the data we interpolate, eg sww info. 533 533 534 534 Output: 535 535 Interpolated values at inputted points. … … 537 537 538 538 interp = Interpolate(vertices, 539 triangles, 539 triangles, 540 540 max_vertices_per_cell=max_points_per_cell, 541 541 mesh_origin=mesh_origin) 542 542 543 543 calc = interp.interpolate(vertex_attributes, 544 544 points, … … 554 554 # @param velocity_y_file Name of the output y velocity file. 555 555 # @param stage_file Name of the output stage file. 556 # @param froude_file 556 # @param froude_file 557 557 # @param time_thinning Time thinning step to use. 558 558 # @param verbose True if this function is to be verbose. … … 604 604 time_thinning=time_thinning, 605 605 use_cache=use_cache) 606 606 607 607 depth_writer = writer(file(depth_file, "wb")) 608 608 velocity_x_writer = writer(file(velocity_x_file, "wb")) … … 620 620 velocity_y_writer.writerow(heading) 621 621 if stage_file is not None: 622 stage_writer.writerow(heading) 622 stage_writer.writerow(heading) 623 623 if froude_file is not None: 624 froude_writer.writerow(heading) 625 624 froude_writer.writerow(heading) 625 626 626 for time in callable_sww.get_time(): 627 627 depths = [time] 628 628 velocity_xs = [time] 629 629 velocity_ys = [time] 630 if stage_file is not None: 631 stages = [time] 632 if froude_file is not None: 633 froudes = [time] 630 if stage_file is not None: 631 stages = [time] 632 if froude_file is not None: 633 froudes = [time] 634 634 for point_i, point in enumerate(points): 635 635 quantities = callable_sww(time,point_i) 636 636 637 637 w = quantities[0] 638 638 z = quantities[1] 639 639 momentum_x = quantities[2] 640 640 momentum_y = quantities[3] 641 depth = w - z 642 641 depth = w - z 642 643 643 if w == NAN or z == NAN or momentum_x == NAN: 644 644 velocity_x = NAN … … 677 677 678 678 if stage_file is not None: 679 stage_writer.writerow(stages) 679 stage_writer.writerow(stages) 680 680 if froude_file is not None: 681 froude_writer.writerow(froudes) 681 froude_writer.writerow(froudes) 682 682 683 683 684 684 ## 685 # @brief 685 # @brief 686 686 class Interpolation_function: 687 687 """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) … … 695 695 Mandatory input 696 696 time: px1 array of monotonously increasing times (float) 697 quantities: Dictionary of arrays or 1 array (float) 697 quantities: Dictionary of arrays or 1 array (float) 698 698 The arrays must either have dimensions pxm or mx1. 699 699 The resulting function will be time dependent in 700 700 the former case while it will be constant with 701 701 respect to time in the latter case. 702 702 703 703 Optional input: 704 704 quantity_names: List of keys into the quantities dictionary for … … 706 706 vertex_coordinates: mx2 array of coordinates (float) 707 707 triangles: nx3 array of indices into vertex_coordinates (Int) 708 interpolation_points: Nx2 array of coordinates to be interpolated to 708 interpolation_points: Nx2 array of coordinates to be interpolated to 709 709 verbose: Level of reporting 710 710 … … 742 742 time, 743 743 quantities, 744 quantity_names=None, 744 quantity_names=None, 745 745 vertex_coordinates=None, 746 746 triangles=None, … … 762 762 print 'Interpolation_function: input checks' 763 763 764 764 # Check temporal info 765 765 time = ensure_numeric(time) 766 766 if not num.alltrue(time[1:] - time[:-1] >= 0): … … 791 791 if triangles is not None: 792 792 triangles = ensure_numeric(triangles) 793 self.spatial = True 793 self.spatial = True 794 794 795 795 if verbose is True: 796 796 print 'Interpolation_function: thinning by %d' % time_thinning 797 797 798 798 799 799 # Thin timesteps if needed 800 800 # Note array() is used to make the thinned arrays contiguous in memory 801 self.time = num.array(time[::time_thinning]) 801 self.time = num.array(time[::time_thinning]) 802 802 for name in quantity_names: 803 803 if len(quantities[name].shape) == 2: … … 806 806 if verbose is True: 807 807 print 'Interpolation_function: precomputing' 808 808 809 809 # Save for use with statistics 810 810 self.quantities_range = {} … … 812 812 q = quantities[name][:].flatten() 813 813 self.quantities_range[name] = [min(q), max(q)] 814 815 self.quantity_names = quantity_names 816 self.vertex_coordinates = vertex_coordinates 814 815 self.quantity_names = quantity_names 816 self.vertex_coordinates = vertex_coordinates 817 817 self.interpolation_points = interpolation_points 818 818 … … 825 825 #if self.spatial is False: 826 826 # raise 'Triangles and vertex_coordinates must be specified' 827 # 827 # 828 828 try: 829 829 self.interpolation_points = \ 830 830 interpolation_points = ensure_numeric(interpolation_points) 831 831 except: 832 832 msg = 'Interpolation points must be an N x 2 numeric array ' \ 833 833 'or a list of points\n' 834 834 msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...') 835 835 raise msg 836 836 … … 839 839 # mesh boundary as defined by triangles and vertex_coordinates. 840 840 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 841 from anuga.utilities.polygon import outside_polygon 841 from anuga.utilities.polygon import outside_polygon 842 842 843 843 # Create temporary mesh object from mesh info passed 844 # into this function. 844 # into this function. 845 845 mesh = Mesh(vertex_coordinates, triangles) 846 846 mesh_boundary_polygon = mesh.get_boundary_polygon() … … 868 868 # FIXME (Ole): Why only Windoze? 869 869 from anuga.utilities.polygon import plot_polygons 870 #out_interp_pts = num.take(interpolation_points,[indices])871 870 title = ('Interpolation points fall ' 872 871 'outside specified mesh') … … 886 885 print msg 887 886 #raise Exception(msg) 888 887 889 888 elif triangles is None and vertex_coordinates is not None: #jj 890 889 #Dealing with sts file … … 911 910 m = len(self.interpolation_points) 912 911 p = len(self.time) 913 914 912 913 for name in quantity_names: 915 914 self.precomputed_values[name] = num.zeros((p, m), num.float) 916 915 … … 918 917 print 'Build interpolator' 919 918 920 919 921 920 # Build interpolator 922 921 if triangles is not None and vertex_coordinates is not None: 923 if verbose: 922 if verbose: 924 923 msg = 'Building interpolation matrix from source mesh ' 925 924 msg += '(%d vertices, %d triangles)' \ … … 927 926 triangles.shape[0]) 928 927 print msg 929 928 930 929 # This one is no longer needed for STS files 931 930 interpol = Interpolate(vertex_coordinates, 932 931 triangles, 933 verbose=verbose) 934 932 verbose=verbose) 933 935 934 elif triangles is None and vertex_coordinates is not None: 936 935 if verbose: … … 943 942 print 'Interpolating (%d interpolation points, %d timesteps).' \ 944 943 %(self.interpolation_points.shape[0], self.time.shape[0]), 945 944 946 945 if time_thinning > 1: 947 946 print 'Timesteps were thinned by a factor of %d' \ … … 950 949 print 951 950 952 951 for i, t in enumerate(self.time): 953 952 # Interpolate quantities at this timestep 954 953 if verbose and i%((p+10)/10) == 0: 955 954 print ' time step %d of %d' %(i, p) 956 955 957 956 for name in quantity_names: 958 957 if len(quantities[name].shape) == 2: … … 963 962 if verbose and i%((p+10)/10) == 0: 964 963 print ' quantity %s, size=%d' %(name, len(Q)) 965 966 # Interpolate 964 965 # Interpolate 967 966 if triangles is not None and vertex_coordinates is not None: 968 967 result = interpol.interpolate(Q, … … 976 975 interpolation_points=\ 977 976 self.interpolation_points) 978 977 979 978 #assert len(result), len(interpolation_points) 980 979 self.precomputed_values[name][i, :] = result … … 1003 1002 def __call__(self, t, point_id=None, x=None, y=None): 1004 1003 """Evaluate f(t) or f(t, point_id) 1005 1006 1007 1008 1009 1010 1011 are None an exception is raised 1012 1004 1005 Inputs: 1006 t: time - Model time. Must lie within existing timesteps 1007 point_id: index of one of the preprocessed points. 1008 1009 If spatial info is present and all of point_id 1010 are None an exception is raised 1011 1013 1012 If no spatial info is present, point_id arguments are ignored 1014 1013 making f a function of time only. 1015 1014 1016 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1017 1018 1015 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1016 FIXME: point_id could also be a slice 1017 FIXME: What if x and y are vectors? 1019 1018 FIXME: What about f(x,y) without t? 1020 1019 """ 1021 1020 1022 1021 from math import pi, cos, sin, sqrt 1023 from anuga.abstract_2d_finite_volumes.util import mean 1022 from anuga.abstract_2d_finite_volumes.util import mean 1024 1023 1025 1024 if self.spatial is True: … … 1057 1056 # Compute interpolated values 1058 1057 q = num.zeros(len(self.quantity_names), num.float) 1059 1058 for i, name in enumerate(self.quantity_names): 1060 1059 Q = self.precomputed_values[name] 1061 1060 1062 1061 if self.spatial is False: 1063 # If there is no spatial info 1062 # If there is no spatial info 1064 1063 assert len(Q.shape) == 1 1065 1064 … … 1076 1075 Q1 = Q[self.index+1, point_id] 1077 1076 1078 # Linear temporal interpolation 1077 # Linear temporal interpolation 1079 1078 if ratio > 0: 1080 1079 if Q0 == NAN and Q1 == NAN: … … 1092 1091 # Replicate q according to x and y 1093 1092 # This is e.g used for Wind_stress 1094 if x is None or y is None: 1093 if x is None or y is None: 1095 1094 return q 1096 1095 else: … … 1106 1105 for col in q: 1107 1106 res.append(col*num.ones(N, num.float)) 1108 1107 1109 1108 return res 1110 1109 … … 1122 1121 """Output statistics about interpolation_function 1123 1122 """ 1124 1123 1125 1124 vertex_coordinates = self.vertex_coordinates 1126 interpolation_points = self.interpolation_points 1125 interpolation_points = self.interpolation_points 1127 1126 quantity_names = self.quantity_names 1128 1127 #quantities = self.quantities 1129 precomputed_values = self.precomputed_values 1130 1128 precomputed_values = self.precomputed_values 1129 1131 1130 x = vertex_coordinates[:,0] 1132 y = vertex_coordinates[:,1] 1131 y = vertex_coordinates[:,1] 1133 1132 1134 1133 str = '------------------------------------------------\n' … … 1144 1143 for name in quantity_names: 1145 1144 minq, maxq = self.quantities_range[name] 1146 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1145 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1147 1146 #q = quantities[name][:].flatten() 1148 1147 #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) 1149 1148 1150 if interpolation_points is not None: 1149 if interpolation_points is not None: 1151 1150 str += ' Interpolation points (xi, eta):'\ 1152 1151 ' number of points == %d\n' %interpolation_points.shape[0] … … 1156 1155 max(interpolation_points[:,1])) 1157 1156 str += ' Interpolated quantities (over all timesteps):\n' 1158 1157 1159 1158 for name in quantity_names: 1160 1159 q = precomputed_values[name][:].flatten() … … 1185 1184 print "x",x 1186 1185 print "y",y 1187 1186 1188 1187 print "time", time 1189 1188 print "quantities", quantities … … 1195 1194 interp = Interpolation_interface(time, 1196 1195 quantities, 1197 quantity_names=quantity_names, 1196 quantity_names=quantity_names, 1198 1197 vertex_coordinates=vertex_coordinates, 1199 1198 triangles=volumes, … … 1208 1207 """ 1209 1208 obsolete - Nothing should be calling this 1210 1209 1211 1210 Read in an sww file. 1212 1211 1213 1212 Input; 1214 1213 file_name - the sww file 1215 1214 1216 1215 Output; 1217 1216 x - Vector of x values … … 1226 1225 msg = 'Function read_sww in interpolat.py is obsolete' 1227 1226 raise Exception, msg 1228 1227 1229 1228 #FIXME Have this reader as part of data_manager? 1230 1231 from Scientific.IO.NetCDF import NetCDFFile 1229 1230 from Scientific.IO.NetCDF import NetCDFFile 1232 1231 import tempfile 1233 1232 import sys 1234 1233 import os 1235 1234 1236 1235 #Check contents 1237 1236 #Get NetCDF 1238 1237 1239 1238 # see if the file is there. Throw a QUIET IO error if it isn't 1240 1239 # I don't think I could implement the above 1241 1240 1242 1241 #throws prints to screen if file not present 1243 1242 junk = tempfile.mktemp(".txt") … … 1245 1244 stdout = sys.stdout 1246 1245 sys.stdout = fd 1247 fid = NetCDFFile(file_name, netcdf_mode_r) 1246 fid = NetCDFFile(file_name, netcdf_mode_r) 1248 1247 sys.stdout = stdout 1249 1248 fd.close() 1250 1249 #clean up 1251 os.remove(junk) 1252 1250 os.remove(junk) 1251 1253 1252 # Get the variables 1254 1253 x = fid.variables['x'][:] 1255 1254 y = fid.variables['y'][:] 1256 volumes = fid.variables['volumes'][:] 1255 volumes = fid.variables['volumes'][:] 1257 1256 time = fid.variables['time'][:] 1258 1257 … … 1266 1265 for name in keys: 1267 1266 quantities[name] = fid.variables[name][:] 1268 1267 1269 1268 fid.close() 1270 1269 return x, y, volumes, time, quantities -
branches/numpy/anuga/geospatial_data/geospatial_data.py
r6360 r6410 216 216 else: 217 217 self.data_points = ensure_numeric(data_points) 218 return219 220 print 'self.data_points=%s' % str(self.data_points)221 print 'self.data_points.shape=%s' % str(self.data_points.shape)222 218 if not (0,) == self.data_points.shape: 223 219 assert len(self.data_points.shape) == 2 … … 551 547 552 548 attributes = {} 553 if file_name[-4:] == ".pts":549 if file_name[-4:] == ".pts": 554 550 try: 555 551 data_points, attributes, geo_reference = \ … … 558 554 msg = 'Could not open file %s ' % file_name 559 555 raise IOError, msg 560 elif file_name[-4:] == ".txt" or file_name[-4:]== ".csv":556 elif file_name[-4:] == ".txt" or file_name[-4:]== ".csv": 561 557 try: 562 558 data_points, attributes, geo_reference = \ … … 588 584 def export_points_file(self, file_name, absolute=True, 589 585 as_lat_long=False, isSouthHemisphere=True): 590 591 586 """write a points file as a text (.csv) or binary (.pts) file 592 587 … … 659 654 # FIXME: add the geo_reference to this 660 655 points = self.get_data_points() 661 sampled_points = num.take(points, indices )656 sampled_points = num.take(points, indices, axis=0) 662 657 663 658 attributes = self.get_all_attributes() … … 666 661 if attributes is not None: 667 662 for key, att in attributes.items(): 668 sampled_attributes[key] = num.take(att, indices )663 sampled_attributes[key] = num.take(att, indices, axis=0) 669 664 670 665 return Geospatial_data(sampled_points, sampled_attributes) … … 695 690 """ 696 691 697 i =0692 i = 0 698 693 self_size = len(self) 699 694 random_list = [] … … 704 699 if verbose: print "make unique random number list and get indices" 705 700 706 total =num.array(range(self_size))701 total = num.array(range(self_size), num.int) #array default# 707 702 total_list = total.tolist() 708 703 … … 734 729 if verbose: print "make random number list and get indices" 735 730 736 j =0737 k =1731 j = 0 732 k = 1 738 733 remainder_list = total_list[:] 739 734 … … 1077 1072 """ 1078 1073 1079 line = file_pointer.readline() .strip()1074 line = file_pointer.readline() 1080 1075 header = clean_line(line, delimiter) 1081 1076 … … 1252 1247 # Create new file 1253 1248 outfile.institution = 'Geoscience Australia' 1254 outfile.description = 'NetCDF format for compact and portable storage ' \1255 'of spatial point data'1249 outfile.description = ('NetCDF format for compact and portable storage ' 1250 'of spatial point data') 1256 1251 1257 1252 # Dimension definitions … … 1578 1573 1579 1574 1580 attribute_smoothed ='elevation'1575 attribute_smoothed = 'elevation' 1581 1576 1582 1577 if mesh_file is None: … … 1637 1632 data = num.array([], dtype=num.float) 1638 1633 1639 data =num.resize(data, (len(points), 3+len(alphas)))1634 data = num.resize(data, (len(points), 3+len(alphas))) 1640 1635 1641 1636 # gets relative point from sample … … 1645 1640 data[:,2] = elevation_sample 1646 1641 1647 normal_cov =num.array(num.zeros([len(alphas), 2]), dtype=num.float)1642 normal_cov = num.array(num.zeros([len(alphas), 2]), dtype=num.float) 1648 1643 1649 1644 if verbose: print 'Setup computational domains with different alphas' … … 1684 1679 1685 1680 normal_cov0 = normal_cov[:,0] 1686 normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0) )1681 normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0) 1687 1682 1688 1683 if plot_name is not None: … … 1903 1898 1904 1899 normal_cov0 = normal_cov[:,0] 1905 normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0) )1900 normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0) 1906 1901 1907 1902 if plot_name is not None: -
branches/numpy/anuga/geospatial_data/test_geospatial_data.py
r6360 r6410 207 207 new_geo_ref = Geo_reference(56, x_p, y_p) 208 208 G.set_geo_reference(new_geo_ref) 209 msg = ('points_ab= %s, but\nG.get_data_points(absolute=True)=%s'209 msg = ('points_ab=\n%s\nbut G.get_data_points(absolute=True)=\n%s' 210 210 % (str(points_ab), str(G.get_data_points(absolute=True)))) 211 211 assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg … … 342 342 343 343 assert num.allclose(P, num.concatenate( (P1,P2) )) 344 msg = 'P=%s' % str(P) 344 msg = ('P=\n%s\nshould be close to\n' 345 '[[2.0, 4.1], [4.0, 7.3],\n' 346 ' [5.1, 9.1], [6.1, 6.3]]' 347 % str(P)) 345 348 assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3], 346 349 [5.1, 9.1], [6.1, 6.3]]), msg … … 369 372 370 373 P1 = G1.get_data_points(absolute=False) 371 msg = ('P1= %s, but\npoints1 - <...>=%s'374 msg = ('P1=\n%s\nbut points1 - <...>=\n%s' 372 375 % (str(P1), 373 376 str(points1 - [geo_ref1.get_xllcorner(), … … 439 442 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 440 443 P = G.get_data_points(absolute=True) 441 msg = 'P= %s' % str(P)444 msg = 'P=\n%s' % str(P) 442 445 assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg 443 446 … … 1022 1025 os.remove(fileName) 1023 1026 1024 msg = ('results.get_data_points()= %s, but\nG.get_data_points(True)=%s'1027 msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s' 1025 1028 % (str(results.get_data_points()), str(G.get_data_points(True)))) 1026 1029 assert num.allclose(results.get_data_points(), -
branches/numpy/anuga/load_mesh/loadASCII.py
r6360 r6410 626 626 num.array(mesh['vertex_attribute_titles'], num.character) 627 627 mesh['segments'] = num.array(mesh['segments'], IntType) 628 print ("Before num.array(), mesh['segment_tags']=%s"629 % str(mesh['segment_tags']))630 628 mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character) 631 print ("After num.array(), mesh['segment_tags']=%s"632 % str(mesh['segment_tags']))633 print "mesh['segment_tags'].shape=%s" % str(mesh['segment_tags'].shape)634 629 mesh['triangles'] = num.array(mesh['triangles'], IntType) 635 630 mesh['triangle_tags'] = num.array(mesh['triangle_tags']) … … 1106 1101 def take_points(dict, indices_to_keep): 1107 1102 dict = point_atts2array(dict) 1108 dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep )1103 dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0) 1109 1104 1110 1105 for key in dict['attributelist'].keys(): 1111 1106 dict['attributelist'][key] = num.take(dict['attributelist'][key], 1112 indices_to_keep )1107 indices_to_keep, axis=0) 1113 1108 1114 1109 return dict -
branches/numpy/anuga/load_mesh/test_loadASCII.py
r6360 r6410 472 472 if __name__ == '__main__': 473 473 suite = unittest.makeSuite(loadASCIITestCase,'test') 474 #suite = unittest.makeSuite(loadASCIITestCase,'test_read_write_msh_file6')475 474 runner = unittest.TextTestRunner() #verbosity=2) 476 475 runner.run(suite) -
branches/numpy/anuga/mesh_engine/mesh_engine_c_layer.c
r6304 r6410 60 60 #include "Python.h" 61 61 62 #include "util_ext.h" 63 64 62 65 //#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH 63 66 … … 119 122 return NULL; 120 123 } 124 125 // check that numpy array objects arrays are C contiguous memory 126 CHECK_C_CONTIG(pointlist); 127 CHECK_C_CONTIG(seglist); 128 CHECK_C_CONTIG(holelist); 129 CHECK_C_CONTIG(regionlist); 130 CHECK_C_CONTIG(pointattributelist); 131 CHECK_C_CONTIG(segmarkerlist); 121 132 122 133 /* Initialize points */ -
branches/numpy/anuga/mesh_engine/test_all.py
r4898 r6410 76 76 return unittest.TestSuite(map(load, modules)) 77 77 78 78 79 if __name__ == '__main__': 79 80 suite = regressionTest() -
branches/numpy/anuga/mesh_engine/test_generate_mesh.py
r6304 r6410 463 463 correct.flat), 464 464 'Failed') 465 465 466 466 467 if __name__ == "__main__": 467 468 468 suite = unittest.makeSuite(triangTestCase,'test') 469 469 runner = unittest.TextTestRunner() #verbosity=2) -
branches/numpy/anuga/pmesh/test_mesh.py
r6360 r6410 2316 2316 if __name__ == "__main__": 2317 2317 suite = unittest.makeSuite(meshTestCase,'test') 2318 #suite = unittest.makeSuite(meshTestCase,'test_import_mesh')2319 2318 runner = unittest.TextTestRunner() #verbosity=2) 2320 2319 runner.run(suite) -
branches/numpy/anuga/shallow_water/data_manager.py
r6304 r6410 2669 2669 # Interpolate using quantity values 2670 2670 if verbose: print 'Interpolating' 2671 interpolated_values = interp.interpolate(q, data_points).flatten 2671 interpolated_values = interp.interpolate(q, data_points).flatten() #????# 2672 2672 2673 2673 if verbose: … … 7308 7308 # FIXME (Ole): Make a generic polygon input check in polygon.py 7309 7309 # and call it here 7310 7311 points = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), axis=1)7312 7310 points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis], 7311 y[:,num.newaxis]), 7312 axis=1)) 7313 7313 point_indices = inside_polygon(points, polygon) 7314 7314 7315 7315 # Restrict quantities to polygon 7316 elevation = num.take(elevation, point_indices )7316 elevation = num.take(elevation, point_indices, axis=0) 7317 7317 stage = num.take(stage, point_indices, axis=1) 7318 7318 7319 7319 # Get info for location of maximal runup 7320 points_in_polygon = num.take(points, point_indices) 7320 points_in_polygon = num.take(points, point_indices, axis=0) 7321 7321 7322 x = points_in_polygon[:,0] 7322 7323 y = points_in_polygon[:,1] … … 7376 7377 else: 7377 7378 # Find maximum elevation among wet nodes 7378 wet_elevation = num.take(elevation, wet_nodes )7379 wet_elevation = num.take(elevation, wet_nodes, axis=0) 7379 7380 runup_index = num.argmax(wet_elevation) 7380 7381 runup = max(wet_elevation) … … 7385 7386 7386 7387 # Record location 7387 wet_x = num.take(x, wet_nodes )7388 wet_y = num.take(y, wet_nodes )7388 wet_x = num.take(x, wet_nodes, axis=0) 7389 wet_y = num.take(y, wet_nodes, axis=0) 7389 7390 maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]] 7390 7391 -
branches/numpy/anuga/shallow_water/shallow_water_domain.py
r6304 r6410 1 """Finite-volume computations of the shallow water wave equation. 1 """ 2 Finite-volume computations of the shallow water wave equation. 2 3 3 4 Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume … … 53 54 Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999 54 55 55 Hydrodynamic modelling of coastal inundation. 56 Hydrodynamic modelling of coastal inundation. 56 57 Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman 57 58 In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on … … 71 72 $Author$ 72 73 $Date$ 73 74 74 """ 75 75 … … 79 79 # $LastChangedRevision$ 80 80 # $LastChangedBy$ 81 81 82 82 83 import numpy as num … … 109 110 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 110 111 111 from anuga.fit_interpolate.interpolate import Modeltime_too_late, Modeltime_too_early 112 113 from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon 112 from anuga.fit_interpolate.interpolate import Modeltime_too_late, \ 113 Modeltime_too_early 114 115 from anuga.utilities.polygon import inside_polygon, polygon_area, \ 116 is_inside_polygon 114 117 115 118 from types import IntType, FloatType … … 117 120 118 121 119 120 # 122 ################################################################################ 121 123 # Shallow water domain 122 #--------------------- 124 ################################################################################ 125 126 ## 127 # @brief Class for a shallow water domain. 123 128 class Domain(Generic_Domain): 124 129 125 130 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 126 131 other_quantities = ['elevation', 'friction'] 127 132 133 ## 134 # @brief Instantiate a shallow water domain. 135 # @param coordinates 136 # @param vertices 137 # @param boundary 138 # @param tagged_elements 139 # @param geo_reference 140 # @param use_inscribed_circle 141 # @param mesh_filename 142 # @param use_cache 143 # @param verbose 144 # @param full_send_dict 145 # @param ghost_recv_dict 146 # @param processor 147 # @param numproc 148 # @param number_of_full_nodes 149 # @param number_of_full_triangles 128 150 def __init__(self, 129 151 coordinates=None, … … 142 164 number_of_full_nodes=None, 143 165 number_of_full_triangles=None): 144 145 166 146 167 other_quantities = ['elevation', 'friction'] … … 162 183 numproc, 163 184 number_of_full_nodes=number_of_full_nodes, 164 number_of_full_triangles=number_of_full_triangles) 165 185 number_of_full_triangles=number_of_full_triangles) 166 186 167 187 self.set_minimum_allowed_height(minimum_allowed_height) 168 188 169 189 self.maximum_allowed_speed = maximum_allowed_speed 170 190 self.g = g 171 self.beta_w 172 self.beta_w_dry 173 self.beta_uh 191 self.beta_w = beta_w 192 self.beta_w_dry = beta_w_dry 193 self.beta_uh = beta_uh 174 194 self.beta_uh_dry = beta_uh_dry 175 self.beta_vh 195 self.beta_vh = beta_vh 176 196 self.beta_vh_dry = beta_vh_dry 177 197 self.alpha_balance = alpha_balance … … 188 208 self.set_store_vertices_uniquely(False) 189 209 self.minimum_storable_height = minimum_storable_height 190 self.quantities_to_be_stored = ['stage', 'xmomentum','ymomentum']210 self.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] 191 211 192 212 # Limiters … … 195 215 self.use_centroid_velocities = use_centroid_velocities 196 216 197 217 ## 218 # @brief 219 # @param beta 198 220 def set_all_limiters(self, beta): 199 """Shorthand to assign one constant value [0,1 [to all limiters.221 """Shorthand to assign one constant value [0,1] to all limiters. 200 222 0 Corresponds to first order, where as larger values make use of 201 the second order scheme. 202 """ 203 204 self.beta_w 205 self.beta_w_dry 223 the second order scheme. 224 """ 225 226 self.beta_w = beta 227 self.beta_w_dry = beta 206 228 self.quantities['stage'].beta = beta 207 208 self.beta_uh 229 230 self.beta_uh = beta 209 231 self.beta_uh_dry = beta 210 232 self.quantities['xmomentum'].beta = beta 211 212 self.beta_vh 233 234 self.beta_vh = beta 213 235 self.beta_vh_dry = beta 214 236 self.quantities['ymomentum'].beta = beta 215 216 217 237 238 ## 239 # @brief 240 # @param flag 241 # @param reduction 218 242 def set_store_vertices_uniquely(self, flag, reduction=None): 219 243 """Decide whether vertex values should be stored uniquely as … … 230 254 #self.reduction = min #Looks better near steep slopes 231 255 232 256 ## 257 # @brief Set the minimum depth that will be written to an SWW file. 258 # @param minimum_storable_height The minimum stored height (in m). 233 259 def set_minimum_storable_height(self, minimum_storable_height): 234 """ 235 Set the minimum depth that will be recognised when writing 260 """Set the minimum depth that will be recognised when writing 236 261 to an sww file. This is useful for removing thin water layers 237 262 that seems to be caused by friction creep. … … 239 264 The minimum allowed sww depth is in meters. 240 265 """ 266 241 267 self.minimum_storable_height = minimum_storable_height 242 268 243 269 ## 270 # @brief 271 # @param minimum_allowed_height 244 272 def set_minimum_allowed_height(self, minimum_allowed_height): 245 """ 246 Set the minimum depth that will be recognised in the numerical 247 scheme 273 """Set the minimum depth recognised in the numerical scheme. 248 274 249 275 The minimum allowed depth is in meters. … … 258 284 #and deal with them adaptively similarly to how we used to use 1 order 259 285 #steps to recover. 286 260 287 self.minimum_allowed_height = minimum_allowed_height 261 self.H0 = minimum_allowed_height 262 263 288 self.H0 = minimum_allowed_height 289 290 ## 291 # @brief 292 # @param maximum_allowed_speed 264 293 def set_maximum_allowed_speed(self, maximum_allowed_speed): 265 """ 266 Set the maximum particle speed that is allowed in water 294 """Set the maximum particle speed that is allowed in water 267 295 shallower than minimum_allowed_height. This is useful for 268 296 controlling speeds in very thin layers of water and at the same time 269 297 allow some movement avoiding pooling of water. 270 271 """ 298 """ 299 272 300 self.maximum_allowed_speed = maximum_allowed_speed 273 301 274 275 def set_points_file_block_line_size(self,points_file_block_line_size): 276 """ 277 Set the minimum depth that will be recognised when writing 302 ## 303 # @brief 304 # @param points_file_block_line_size 305 def set_points_file_block_line_size(self, points_file_block_line_size): 306 """Set the minimum depth that will be recognised when writing 278 307 to an sww file. This is useful for removing thin water layers 279 308 that seems to be caused by friction creep. … … 282 311 """ 283 312 self.points_file_block_line_size = points_file_block_line_size 284 285 313 314 ## 315 # @brief Set the quantities that will be written to an SWW file. 316 # @param q The quantities to be written. 317 # @note Param 'q' may be None, single quantity or list of quantity strings. 318 # @note If 'q' is None, no quantities will be stored in the SWW file. 286 319 def set_quantities_to_be_stored(self, q): 287 320 """Specify which quantities will be stored in the sww file. … … 295 328 each yieldstep (This is in addition to the quantities elevation 296 329 and friction) 297 330 298 331 If q is None, storage will be switched off altogether. 299 332 """ 300 301 333 302 334 if q is None: … … 310 342 # Check correcness 311 343 for quantity_name in q: 312 msg = 'Quantity %s is not a valid conserved quantity'\ 313 %quantity_name 314 344 msg = ('Quantity %s is not a valid conserved quantity' 345 % quantity_name) 315 346 assert quantity_name in self.conserved_quantities, msg 316 347 317 348 self.quantities_to_be_stored = q 318 349 319 320 350 ## 351 # @brief 352 # @param indices 321 353 def get_wet_elements(self, indices=None): 322 354 """Return indices for elements where h > minimum_allowed_height … … 328 360 indices = get_wet_elements() 329 361 330 Note, centroid values are used for this operation 362 Note, centroid values are used for this operation 331 363 """ 332 364 … … 335 367 from anuga.config import minimum_allowed_height 336 368 337 338 369 elevation = self.get_quantity('elevation').\ 370 get_values(location='centroids', indices=indices) 371 stage = self.get_quantity('stage').\ 339 372 get_values(location='centroids', indices=indices) 340 stage = self.get_quantity('stage').\341 get_values(location='centroids', indices=indices)342 373 depth = stage - elevation 343 374 … … 345 376 wet_indices = num.compress(depth > minimum_allowed_height, 346 377 num.arange(len(depth))) 347 return wet_indices 348 349 378 return wet_indices 379 380 ## 381 # @brief 382 # @param indices 350 383 def get_maximum_inundation_elevation(self, indices=None): 351 384 """Return highest elevation where h > 0 … … 357 390 q = get_maximum_inundation_elevation() 358 391 359 Note, centroid values are used for this operation 392 Note, centroid values are used for this operation 360 393 """ 361 394 362 395 wet_elements = self.get_wet_elements(indices) 363 396 return self.get_quantity('elevation').\ 364 get_maximum_value(indices=wet_elements) 365 366 397 get_maximum_value(indices=wet_elements) 398 399 ## 400 # @brief 401 # @param indices 367 402 def get_maximum_inundation_location(self, indices=None): 368 403 """Return location of highest elevation where h > 0 … … 374 409 q = get_maximum_inundation_location() 375 410 376 Note, centroid values are used for this operation 411 Note, centroid values are used for this operation 377 412 """ 378 413 379 414 wet_elements = self.get_wet_elements(indices) 380 415 return self.get_quantity('elevation').\ 381 get_maximum_location(indices=wet_elements) 382 383 384 385 386 def get_flow_through_cross_section(self, polyline, 387 verbose=False): 388 """Get the total flow through an arbitrary poly line. 389 390 This is a run-time equivalent of the function with same name 416 get_maximum_location(indices=wet_elements) 417 418 ## 419 # @brief Get the total flow through an arbitrary poly line. 420 # @param polyline Representation of desired cross section. 421 # @param verbose True if this method is to be verbose. 422 # @note 'polyline' may contain multiple sections allowing complex shapes. 423 # @note Assume absolute UTM coordinates. 424 def get_flow_through_cross_section(self, polyline, verbose=False): 425 """Get the total flow through an arbitrary poly line. 426 427 This is a run-time equivalent of the function with same name 391 428 in data_manager.py 392 429 393 430 Input: 394 polyline: Representation of desired cross section - it may contain 395 multiple sections allowing for complex shapes. Assume 431 polyline: Representation of desired cross section - it may contain 432 multiple sections allowing for complex shapes. Assume 396 433 absolute UTM coordinates. 397 Format [[x0, y0], [x1, y1], ...] 398 399 Output: 434 Format [[x0, y0], [x1, y1], ...] 435 436 Output: 400 437 Q: Total flow [m^3/s] across given segments. 401 402 403 """ 404 438 """ 439 405 440 # Find all intersections and associated triangles. 406 segments = self.get_intersecting_segments(polyline, 407 use_cache=True, 441 segments = self.get_intersecting_segments(polyline, use_cache=True, 408 442 verbose=verbose) 409 443 410 444 # Get midpoints 411 midpoints = segment_midpoints(segments) 412 445 midpoints = segment_midpoints(segments) 446 413 447 # Make midpoints Geospatial instances 414 midpoints = ensure_geospatial(midpoints, self.geo_reference) 415 416 # Compute flow 448 midpoints = ensure_geospatial(midpoints, self.geo_reference) 449 450 # Compute flow 417 451 if verbose: print 'Computing flow through specified cross section' 418 452 419 453 # Get interpolated values 420 454 xmomentum = self.get_quantity('xmomentum') 421 ymomentum = self.get_quantity('ymomentum') 422 423 uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True) 424 vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True) 425 455 ymomentum = self.get_quantity('ymomentum') 456 457 uh = xmomentum.get_values(interpolation_points=midpoints, 458 use_cache=True) 459 vh = ymomentum.get_values(interpolation_points=midpoints, 460 use_cache=True) 461 426 462 # Compute and sum flows across each segment 427 total_flow =0463 total_flow = 0 428 464 for i in range(len(uh)): 429 430 # Inner product of momentum vector with segment normal [m^2/s] 465 # Inner product of momentum vector with segment normal [m^2/s] 431 466 normal = segments[i].normal 432 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 433 467 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 468 434 469 # Flow across this segment [m^3/s] 435 470 segment_flow = normal_momentum*segments[i].length … … 437 472 # Accumulate 438 473 total_flow += segment_flow 439 474 440 475 return total_flow 441 476 442 443 477 ## 478 # @brief 479 # @param polyline Representation of desired cross section. 480 # @param kind Select energy type to compute ('specific' or 'total'). 481 # @param verbose True if this method is to be verbose. 482 # @note 'polyline' may contain multiple sections allowing complex shapes. 483 # @note Assume absolute UTM coordinates. 444 484 def get_energy_through_cross_section(self, polyline, 445 485 kind='total', 446 verbose=False): 486 verbose=False): 447 487 """Obtain average energy head [m] across specified cross section. 448 488 449 489 Inputs: 450 polyline: Representation of desired cross section - it may contain 451 multiple sections allowing for complex shapes. Assume 490 polyline: Representation of desired cross section - it may contain 491 multiple sections allowing for complex shapes. Assume 452 492 absolute UTM coordinates. 453 493 Format [[x0, y0], [x1, y1], ...] 454 kind: Select which energy to compute. 494 kind: Select which energy to compute. 455 495 Options are 'specific' and 'total' (default) 456 496 … … 458 498 E: Average energy [m] across given segments for all stored times. 459 499 460 The average velocity is computed for each triangle intersected by 461 the polyline and averaged weighted by segment lengths. 462 463 The typical usage of this function would be to get average energy of 464 flow in a channel, and the polyline would then be a cross section 500 The average velocity is computed for each triangle intersected by 501 the polyline and averaged weighted by segment lengths. 502 503 The typical usage of this function would be to get average energy of 504 flow in a channel, and the polyline would then be a cross section 465 505 perpendicular to the flow. 466 506 467 #FIXME (Ole) - need name for this energy reflecting that its dimension 507 #FIXME (Ole) - need name for this energy reflecting that its dimension 468 508 is [m]. 469 509 """ 470 510 471 from anuga.config import g, epsilon, velocity_protection as h0 472 473 511 from anuga.config import g, epsilon, velocity_protection as h0 512 474 513 # Find all intersections and associated triangles. 475 segments = self.get_intersecting_segments(polyline, 476 use_cache=True, 514 segments = self.get_intersecting_segments(polyline, use_cache=True, 477 515 verbose=verbose) 478 516 479 517 # Get midpoints 480 midpoints = segment_midpoints(segments) 481 518 midpoints = segment_midpoints(segments) 519 482 520 # Make midpoints Geospatial instances 483 midpoints = ensure_geospatial(midpoints, self.geo_reference) 484 485 # Compute energy 486 if verbose: print 'Computing %s energy' %kind 487 521 midpoints = ensure_geospatial(midpoints, self.geo_reference) 522 523 # Compute energy 524 if verbose: print 'Computing %s energy' %kind 525 488 526 # Get interpolated values 489 stage = self.get_quantity('stage') 490 elevation = self.get_quantity('elevation') 527 stage = self.get_quantity('stage') 528 elevation = self.get_quantity('elevation') 491 529 xmomentum = self.get_quantity('xmomentum') 492 ymomentum = self.get_quantity('ymomentum') 530 ymomentum = self.get_quantity('ymomentum') 493 531 494 532 w = stage.get_values(interpolation_points=midpoints, use_cache=True) 495 z = elevation.get_values(interpolation_points=midpoints, use_cache=True) 496 uh = xmomentum.get_values(interpolation_points=midpoints, use_cache=True) 497 vh = ymomentum.get_values(interpolation_points=midpoints, use_cache=True) 498 h = w-z # Depth 499 533 z = elevation.get_values(interpolation_points=midpoints, use_cache=True) 534 uh = xmomentum.get_values(interpolation_points=midpoints, 535 use_cache=True) 536 vh = ymomentum.get_values(interpolation_points=midpoints, 537 use_cache=True) 538 h = w-z # Depth 539 500 540 # Compute total length of polyline for use with weighted averages 501 541 total_line_length = 0.0 502 542 for segment in segments: 503 543 total_line_length += segment.length 504 544 505 545 # Compute and sum flows across each segment 506 average_energy =0.0546 average_energy = 0.0 507 547 for i in range(len(w)): 508 509 548 # Average velocity across this segment 510 549 if h[i] > epsilon: … … 514 553 else: 515 554 u = v = 0.0 516 517 speed_squared = u*u + v*v 555 556 speed_squared = u*u + v*v 518 557 kinetic_energy = 0.5*speed_squared/g 519 558 520 559 if kind == 'specific': 521 560 segment_energy = h[i] + kinetic_energy 522 561 elif kind == 'total': 523 segment_energy = w[i] + kinetic_energy 562 segment_energy = w[i] + kinetic_energy 524 563 else: 525 564 msg = 'Energy kind must be either "specific" or "total".' 526 565 msg += ' I got %s' %kind 527 528 566 529 567 # Add to weighted average 530 568 weigth = segments[i].length/total_line_length 531 569 average_energy += segment_energy*weigth 532 533 570 534 571 return average_energy 535 536 537 538 572 573 ## 574 # @brief Run integrity checks on shallow water domain. 539 575 def check_integrity(self): 540 Generic_Domain.check_integrity(self) 576 Generic_Domain.check_integrity(self) # super integrity check 541 577 542 578 #Check that we are solving the shallow water wave equation 543 544 579 msg = 'First conserved quantity must be "stage"' 545 580 assert self.conserved_quantities[0] == 'stage', msg … … 549 584 assert self.conserved_quantities[2] == 'ymomentum', msg 550 585 586 ## 587 # @brief 551 588 def extrapolate_second_order_sw(self): 552 #Call correct module function 553 #(either from this module or C-extension) 589 #Call correct module function (either from this module or C-extension) 554 590 extrapolate_second_order_sw(self) 555 591 592 ## 593 # @brief 556 594 def compute_fluxes(self): 557 #Call correct module function 558 #(either from this module or C-extension) 595 #Call correct module function (either from this module or C-extension) 559 596 compute_fluxes(self) 560 597 598 ## 599 # @brief 561 600 def distribute_to_vertices_and_edges(self): 562 601 # Call correct module function 563 602 if self.use_edge_limiter: 564 distribute_using_edge_limiter(self) 603 distribute_using_edge_limiter(self) 565 604 else: 566 605 distribute_using_vertex_limiter(self) 567 606 568 569 570 607 ## 608 # @brief Evolve the model by one step. 609 # @param yieldstep 610 # @param finaltime 611 # @param duration 612 # @param skip_initial_step 571 613 def evolve(self, 572 yieldstep = None, 573 finaltime = None, 574 duration = None, 575 skip_initial_step = False): 576 """Specialisation of basic evolve method from parent class 577 """ 614 yieldstep=None, 615 finaltime=None, 616 duration=None, 617 skip_initial_step=False): 618 """Specialisation of basic evolve method from parent class""" 578 619 579 620 # Call check integrity here rather than from user scripts 580 621 # self.check_integrity() 581 622 582 msg = ' Parameter beta_w must be in the interval [0, 2['623 msg = 'Attribute self.beta_w must be in the interval [0, 2]' 583 624 assert 0 <= self.beta_w <= 2.0, msg 584 625 585 586 626 # Initial update of vertex and edge values before any STORAGE 587 # and or visualisation 627 # and or visualisation. 588 628 # This is done again in the initialisation of the Generic_Domain 589 # evolve loop but we do it here to ensure the values are ok for storage 629 # evolve loop but we do it here to ensure the values are ok for storage. 590 630 self.distribute_to_vertices_and_edges() 591 631 592 632 if self.store is True and self.time == 0.0: 593 633 self.initialise_storage() 594 # print 'Storing results in ' + self.writer.filename595 634 else: 596 635 pass … … 601 640 602 641 # Call basic machinery from parent class 603 for t in Generic_Domain.evolve(self, 604 yieldstep=yieldstep, 605 finaltime=finaltime, 606 duration=duration, 642 for t in Generic_Domain.evolve(self, yieldstep=yieldstep, 643 finaltime=finaltime, duration=duration, 607 644 skip_initial_step=skip_initial_step): 608 609 645 # Store model data, e.g. for subsequent visualisation 610 646 if self.store is True: … … 617 653 yield(t) 618 654 619 655 ## 656 # @brief 620 657 def initialise_storage(self): 621 658 """Create and initialise self.writer object for storing data. … … 631 668 self.writer.store_connectivity() 632 669 633 670 ## 671 # @brief 672 # @param name 634 673 def store_timestep(self, name): 635 674 """Store named quantity and time. … … 638 677 self.write has been initialised 639 678 """ 679 640 680 self.writer.store_timestep(name) 641 681 642 682 ## 683 # @brief Get time stepping statistics string for printing. 684 # @param track_speeds 685 # @param triangle_id 643 686 def timestepping_statistics(self, 644 687 track_speeds=False, 645 triangle_id=None): 688 triangle_id=None): 646 689 """Return string with time stepping statistics for printing or logging 647 690 … … 651 694 """ 652 695 653 from anuga.config import epsilon, g 654 696 from anuga.config import epsilon, g 655 697 656 698 # Call basic machinery from parent class 657 msg = Generic_Domain.timestepping_statistics(self, 658 track_speeds, 699 msg = Generic_Domain.timestepping_statistics(self, track_speeds, 659 700 triangle_id) 660 701 661 702 if track_speeds is True: 662 663 703 # qwidth determines the text field used for quantities 664 704 qwidth = self.qwidth 665 705 666 706 # Selected triangle 667 707 k = self.k … … 669 709 # Report some derived quantities at vertices, edges and centroid 670 710 # specific to the shallow water wave equation 671 672 711 z = self.quantities['elevation'] 673 w = self.quantities['stage'] 712 w = self.quantities['stage'] 674 713 675 714 Vw = w.get_values(location='vertices', indices=[k])[0] … … 679 718 Vz = z.get_values(location='vertices', indices=[k])[0] 680 719 Ez = z.get_values(location='edges', indices=[k])[0] 681 Cz = z.get_values(location='centroids', indices=[k]) 682 720 Cz = z.get_values(location='centroids', indices=[k]) 683 721 684 722 name = 'depth' … … 686 724 Eh = Ew-Ez 687 725 Ch = Cw-Cz 688 726 689 727 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 690 728 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2]) 691 729 692 730 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 693 731 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2]) 694 732 695 733 s += ' %s: centroid_value = %.4f\n'\ 696 %(name.ljust(qwidth), Ch[0]) 697 734 %(name.ljust(qwidth), Ch[0]) 735 698 736 msg += s 699 737 … … 704 742 Euh = uh.get_values(location='edges', indices=[k])[0] 705 743 Cuh = uh.get_values(location='centroids', indices=[k]) 706 744 707 745 Vvh = vh.get_values(location='vertices', indices=[k])[0] 708 746 Evh = vh.get_values(location='edges', indices=[k])[0] … … 712 750 Vu = Vuh/(Vh + epsilon) 713 751 Eu = Euh/(Eh + epsilon) 714 Cu = Cuh/(Ch + epsilon) 752 Cu = Cuh/(Ch + epsilon) 715 753 name = 'U' 716 754 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 717 755 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2]) 718 756 719 757 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 720 758 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2]) 721 759 722 760 s += ' %s: centroid_value = %.4f\n'\ 723 %(name.ljust(qwidth), Cu[0]) 724 761 %(name.ljust(qwidth), Cu[0]) 762 725 763 msg += s 726 764 727 765 Vv = Vvh/(Vh + epsilon) 728 766 Ev = Evh/(Eh + epsilon) 729 Cv = Cvh/(Ch + epsilon) 767 Cv = Cvh/(Ch + epsilon) 730 768 name = 'V' 731 769 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 732 770 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2]) 733 771 734 772 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 735 773 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2]) 736 774 737 775 s += ' %s: centroid_value = %.4f\n'\ 738 %(name.ljust(qwidth), Cv[0]) 739 776 %(name.ljust(qwidth), Cv[0]) 777 740 778 msg += s 741 742 779 743 780 # Froude number in each direction … … 746 783 Efx = Eu/(num.sqrt(g*Eh) + epsilon) 747 784 Cfx = Cu/(num.sqrt(g*Ch) + epsilon) 748 785 749 786 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 750 787 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2]) 751 788 752 789 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 753 790 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2]) 754 791 755 792 s += ' %s: centroid_value = %.4f\n'\ 756 %(name.ljust(qwidth), Cfx[0]) 757 793 %(name.ljust(qwidth), Cfx[0]) 794 758 795 msg += s 759 760 796 761 797 name = 'Froude (y)' … … 763 799 Efy = Ev/(num.sqrt(g*Eh) + epsilon) 764 800 Cfy = Cv/(num.sqrt(g*Ch) + epsilon) 765 801 766 802 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 767 803 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2]) 768 804 769 805 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 770 806 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2]) 771 807 772 808 s += ' %s: centroid_value = %.4f\n'\ 773 %(name.ljust(qwidth), Cfy[0]) 774 775 msg += s 776 777 809 %(name.ljust(qwidth), Cfy[0]) 810 811 msg += s 778 812 779 813 return msg 780 781 782 783 #=============== End of class Shallow Water Domain =============================== 784 814 815 ################################################################################ 816 # End of class Shallow Water Domain 817 ################################################################################ 785 818 786 819 #----------------- … … 788 821 #----------------- 789 822 823 ## @brief Compute fluxes and timestep suitable for all volumes in domain. 824 # @param domain The domain to calculate fluxes for. 790 825 def compute_fluxes(domain): 791 """Compute all fluxes and the timestep suitable for all volumes 792 in domain. 826 """Compute fluxes and timestep suitable for all volumes in domain. 793 827 794 828 Compute total flux for each conserved quantity using "flux_function" … … 806 840 domain.explicit_update is reset to computed flux values 807 841 domain.timestep is set to the largest step satisfying all volumes. 808 809 842 810 843 This wrapper calls the underlying C version of compute fluxes … … 812 845 813 846 import sys 814 815 N = len(domain) # number_of_triangles 847 from shallow_water_ext import compute_fluxes_ext_central \ 848 as compute_fluxes_ext 849 850 N = len(domain) # number_of_triangles 816 851 817 852 # Shortcuts … … 822 857 823 858 timestep = float(sys.maxint) 824 from shallow_water_ext import\825 compute_fluxes_ext_central as compute_fluxes_ext826 827 859 828 860 flux_timestep = compute_fluxes_ext(timestep, … … 853 885 domain.flux_timestep = flux_timestep 854 886 855 856 857 #--------------------------------------- 887 ################################################################################ 858 888 # Module functions for gradient limiting 859 #--------------------------------------- 860 861 862 # MH090605 The following method belongs to the shallow_water domain class 863 # see comments in the corresponding method in shallow_water_ext.c 889 ################################################################################ 890 891 ## 892 # @brief Wrapper for C version of extrapolate_second_order_sw. 893 # @param domain The domain to operate on. 894 # @note MH090605 The following method belongs to the shallow_water domain class 895 # see comments in the corresponding method in shallow_water_ext.c 864 896 def extrapolate_second_order_sw(domain): 865 """Wrapper calling C version of extrapolate_second_order_sw 866 """ 897 """Wrapper calling C version of extrapolate_second_order_sw""" 898 867 899 import sys 900 from shallow_water_ext import extrapolate_second_order_sw as extrapol2 868 901 869 902 N = len(domain) # number_of_triangles … … 875 908 Elevation = domain.quantities['elevation'] 876 909 877 from shallow_water_ext import extrapolate_second_order_sw as extrapol2878 910 extrapol2(domain, 879 911 domain.surrogate_neighbours, … … 891 923 int(domain.optimise_dry_cells)) 892 924 893 925 ## 926 # @brief Distribution from centroids to vertices specific to the SWW eqn. 927 # @param domain The domain to operate on. 894 928 def distribute_using_vertex_limiter(domain): 895 """Distribution from centroids to vertices specific to the 896 shallow water wave 897 equation. 929 """Distribution from centroids to vertices specific to the SWW equation. 898 930 899 931 It will ensure that h (w-z) is always non-negative even in the … … 912 944 Postcondition 913 945 Conserved quantities defined at vertices 914 915 946 """ 916 917 918 947 919 948 # Remove very thin layers of water … … 946 975 raise 'Unknown order' 947 976 948 949 977 # Take bed elevation into account when water heights are small 950 978 balance_deep_and_shallow(domain) … … 955 983 Q.interpolate_from_vertices_to_edges() 956 984 957 958 985 ## 986 # @brief Distribution from centroids to edges specific to the SWW eqn. 987 # @param domain The domain to operate on. 959 988 def distribute_using_edge_limiter(domain): 960 """Distribution from centroids to edges specific to the 961 shallow water wave 962 equation. 989 """Distribution from centroids to edges specific to the SWW eqn. 963 990 964 991 It will ensure that h (w-z) is always non-negative even in the … … 976 1003 Postcondition 977 1004 Conserved quantities defined at vertices 978 979 1005 """ 980 1006 981 1007 # Remove very thin layers of water 982 1008 protect_against_infinitesimal_and_negative_heights(domain) 983 984 1009 985 1010 for name in domain.conserved_quantities: … … 999 1024 Q.interpolate_from_vertices_to_edges() 1000 1025 1001 1026 ## 1027 # @brief Protect against infinitesimal heights and associated high velocities. 1028 # @param domain The domain to operate on. 1002 1029 def protect_against_infinitesimal_and_negative_heights(domain): 1003 """Protect against infinitesimal heights and associated high velocities 1004 """ 1030 """Protect against infinitesimal heights and associated high velocities""" 1031 1032 from shallow_water_ext import protect 1005 1033 1006 1034 # Shortcuts … … 1010 1038 ymomc = domain.quantities['ymomentum'].centroid_values 1011 1039 1012 from shallow_water_ext import protect1013 1014 1040 protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 1015 1041 domain.epsilon, wc, zc, xmomc, ymomc) 1016 1042 1017 1018 1043 ## 1044 # @brief Wrapper for C function balance_deep_and_shallow_c(). 1045 # @param domain The domain to operate on. 1019 1046 def balance_deep_and_shallow(domain): 1020 1047 """Compute linear combination between stage as computed by … … 1031 1058 """ 1032 1059 1033 from shallow_water_ext import balance_deep_and_shallow as balance_deep_and_shallow_c1034 1060 from shallow_water_ext import balance_deep_and_shallow \ 1061 as balance_deep_and_shallow_c 1035 1062 1036 1063 # Shortcuts 1037 1064 wc = domain.quantities['stage'].centroid_values 1038 1065 zc = domain.quantities['elevation'].centroid_values 1039 1040 1066 wv = domain.quantities['stage'].vertex_values 1041 1067 zv = domain.quantities['elevation'].vertex_values … … 1050 1076 1051 1077 balance_deep_and_shallow_c(domain, 1052 wc, zc, wv, zv, wc, 1078 wc, zc, wv, zv, wc, 1053 1079 xmomc, ymomc, xmomv, ymomv) 1054 1080 1055 1081 1056 1057 1058 #------------------------------------------------------------------ 1082 ################################################################################ 1059 1083 # Boundary conditions - specific to the shallow water wave equation 1060 #------------------------------------------------------------------ 1084 ################################################################################ 1085 1086 ## 1087 # @brief Class for a reflective boundary. 1088 # @note Inherits from Boundary. 1061 1089 class Reflective_boundary(Boundary): 1062 1090 """Reflective boundary returns same conserved quantities as … … 1068 1096 """ 1069 1097 1070 def __init__(self, domain = None): 1098 ## 1099 # @brief Instantiate a Reflective_boundary. 1100 # @param domain 1101 def __init__(self, domain=None): 1071 1102 Boundary.__init__(self) 1072 1103 1073 1104 if domain is None: 1074 1105 msg = 'Domain must be specified for reflective boundary' 1075 raise msg1106 raise Exception, msg 1076 1107 1077 1108 # Handy shorthands 1078 self.stage 1079 self.xmom 1080 self.ymom 1109 self.stage = domain.quantities['stage'].edge_values 1110 self.xmom = domain.quantities['xmomentum'].edge_values 1111 self.ymom = domain.quantities['ymomentum'].edge_values 1081 1112 self.normals = domain.normals 1082 1113 1083 1114 self.conserved_quantities = num.zeros(3, num.float) 1084 1115 1116 ## 1117 # @brief Return a representation of this instance. 1085 1118 def __repr__(self): 1086 1119 return 'Reflective_boundary' 1087 1120 1088 1121 ## 1122 # @brief Calculate reflections (reverse outward momentum). 1123 # @param vol_id 1124 # @param edge_id 1089 1125 def evaluate(self, vol_id, edge_id): 1090 1126 """Reflective boundaries reverses the outward momentum … … 1099 1135 normal = self.normals[vol_id, 2*edge_id:2*edge_id+2] 1100 1136 1101 1102 1137 r = rotate(q, normal, direction = 1) 1103 1138 r[1] = -r[1] … … 1107 1142 1108 1143 1109 1110 1144 ## 1145 # @brief Class for a transmissive boundary. 1146 # @note Inherits from Boundary. 1111 1147 class Transmissive_momentum_set_stage_boundary(Boundary): 1112 1148 """Returns same momentum conserved quantities as … … 1117 1153 Example: 1118 1154 1119 def waveform(t): 1155 def waveform(t): 1120 1156 return sea_level + normalized_amplitude/cosh(t-25)**2 1121 1157 1122 1158 Bts = Transmissive_momentum_set_stage_boundary(domain, waveform) 1123 1124 1159 1125 1160 Underlying domain must be specified when boundary is instantiated 1126 1161 """ 1127 1162 1128 def __init__(self, domain = None, function=None): 1163 ## 1164 # @brief Instantiate a Reflective_boundary. 1165 # @param domain 1166 # @param function 1167 def __init__(self, domain=None, function=None): 1129 1168 Boundary.__init__(self) 1130 1169 1131 1170 if domain is None: 1132 1171 msg = 'Domain must be specified for this type boundary' 1133 raise msg1172 raise Exception, msg 1134 1173 1135 1174 if function is None: 1136 1175 msg = 'Function must be specified for this type boundary' 1137 raise msg1138 1139 self.domain 1176 raise Exception, msg 1177 1178 self.domain = domain 1140 1179 self.function = function 1141 1180 1181 ## 1182 # @brief Return a representation of this instance. 1142 1183 def __repr__(self): 1143 1184 return 'Transmissive_momentum_set_stage_boundary(%s)' %self.domain 1144 1185 1186 ## 1187 # @brief Calculate transmissive results. 1188 # @param vol_id 1189 # @param edge_id 1145 1190 def evaluate(self, vol_id, edge_id): 1146 1191 """Transmissive momentum set stage boundaries return the edge momentum … … 1150 1195 q = self.domain.get_conserved_quantities(vol_id, edge = edge_id) 1151 1196 1152 1153 1197 t = self.domain.get_time() 1154 1198 1155 1199 if hasattr(self.function, 'time'): 1156 # Roll boundary over if time exceeds 1200 # Roll boundary over if time exceeds 1157 1201 while t > self.function.time[-1]: 1158 msg = 'WARNING: domain time %.2f has exceeded' % t1202 msg = 'WARNING: domain time %.2f has exceeded' % t 1159 1203 msg += 'time provided in ' 1160 1204 msg += 'transmissive_momentum_set_stage boundary object.\n' … … 1163 1207 t -= self.function.time[-1] 1164 1208 1165 1166 1209 value = self.function(t) 1167 1210 try: … … 1169 1212 except: 1170 1213 x = float(value[0]) 1171 1214 1172 1215 q[0] = x 1173 1216 return q 1174 1175 1217 1176 1218 # FIXME: Consider this (taken from File_boundary) to allow … … 1183 1225 1184 1226 1185 # Backward compatibility1186 # FIXME(Ole): Deprecate1227 ## 1228 # @brief Deprecated boundary class. 1187 1229 class Transmissive_Momentum_Set_Stage_boundary(Transmissive_momentum_set_stage_boundary): 1188 1230 pass 1189 1231 1190 1191 1232 1233 ## 1234 # @brief A transmissive boundary, momentum set to zero. 1235 # @note Inherits from Bouondary. 1192 1236 class Transmissive_stage_zero_momentum_boundary(Boundary): 1193 """Return same stage as those present in its neighbour volume. Set momentum to zero. 1237 """Return same stage as those present in its neighbour volume. 1238 Set momentum to zero. 1194 1239 1195 1240 Underlying domain must be specified when boundary is instantiated 1196 1241 """ 1197 1242 1243 ## 1244 # @brief Instantiate a Transmissive (zero momentum) boundary. 1245 # @param domain 1198 1246 def __init__(self, domain=None): 1199 1247 Boundary.__init__(self) 1200 1248 1201 1249 if domain is None: 1202 msg = 'Domain must be specified for '1203 msg += 'Transmissive_stage_zero_momentum boundary'1250 msg = ('Domain must be specified for ' 1251 'Transmissive_stage_zero_momentum boundary') 1204 1252 raise Exception, msg 1205 1253 1206 1254 self.domain = domain 1207 1255 1256 ## 1257 # @brief Return a representation of this instance. 1208 1258 def __repr__(self): 1209 return 'Transmissive_stage_zero_momentum_boundary(%s)' %self.domain 1210 1259 return 'Transmissive_stage_zero_momentum_boundary(%s)' % self.domain 1260 1261 ## 1262 # @brief Calculate transmissive (zero momentum) results. 1263 # @param vol_id 1264 # @param edge_id 1211 1265 def evaluate(self, vol_id, edge_id): 1212 1266 """Transmissive boundaries return the edge values … … 1215 1269 1216 1270 q = self.domain.get_conserved_quantities(vol_id, edge=edge_id) 1217 1271 1218 1272 q[1] = q[2] = 0.0 1219 1273 return q 1220 1274 1221 1275 1222 1276 ## 1277 # @brief Class for a Dirichlet discharge boundary. 1278 # @note Inherits from Boundary. 1223 1279 class Dirichlet_discharge_boundary(Boundary): 1224 1280 """ … … 1229 1285 """ 1230 1286 1287 ## 1288 # @brief Instantiate a Dirichlet discharge boundary. 1289 # @param domain 1290 # @param stage0 1291 # @param wh0 1231 1292 def __init__(self, domain=None, stage0=None, wh0=None): 1232 1293 Boundary.__init__(self) 1233 1294 1234 1295 if domain is None: 1235 msg = 'Domain must be specified for this type boundary'1236 raise msg1296 msg = 'Domain must be specified for this type of boundary' 1297 raise Exception, msg 1237 1298 1238 1299 if stage0 is None: 1239 raise 'set stage'1300 raise Exception, 'Stage must be specified for this type of boundary' 1240 1301 1241 1302 if wh0 is None: 1242 1303 wh0 = 0.0 1243 1304 1244 self.domain 1245 self.stage0 1305 self.domain = domain 1306 self.stage0 = stage0 1246 1307 self.wh0 = wh0 1247 1308 1309 ## 1310 # @brief Return a representation of this instance. 1248 1311 def __repr__(self): 1249 return 'Dirichlet_Discharge_boundary(%s)' %self.domain 1250 1312 return 'Dirichlet_Discharge_boundary(%s)' % self.domain 1313 1314 ## 1315 # @brief Calculate Dirichlet discharge boundary results. 1316 # @param vol_id 1317 # @param edge_id 1251 1318 def evaluate(self, vol_id, edge_id): 1252 """Set discharge in the (inward) normal direction 1253 """ 1319 """Set discharge in the (inward) normal direction""" 1254 1320 1255 1321 normal = self.domain.get_normal(vol_id,edge_id) 1256 1322 q = [self.stage0, -self.wh0*normal[0], -self.wh0*normal[1]] 1257 1323 return q 1258 1259 1324 1260 1325 # FIXME: Consider this (taken from File_boundary) to allow … … 1267 1332 1268 1333 1269 1270 # Backward compatibility 1334 # Backward compatibility 1271 1335 # FIXME(Ole): Deprecate 1336 ## 1337 # @brief Deprecated 1272 1338 class Dirichlet_Discharge_boundary(Dirichlet_discharge_boundary): 1273 1339 pass 1274 1275 1276 1277 1278 1340 1341 1342 ## 1343 # @brief A class for a 'field' boundary. 1344 # @note Inherits from Boundary. 1279 1345 class Field_boundary(Boundary): 1280 """Set boundary from given field represented in an sww file containing values 1281 for stage, xmomentum and ymomentum. 1282 Optionally, the user can specify mean_stage to offset the stage provided in the 1283 sww file. 1284 1285 This function is a thin wrapper around the generic File_boundary. The 1346 """Set boundary from given field represented in an sww file containing 1347 values for stage, xmomentum and ymomentum. 1348 1349 Optionally, the user can specify mean_stage to offset the stage provided 1350 in the sww file. 1351 1352 This function is a thin wrapper around the generic File_boundary. The 1286 1353 difference between the file_boundary and field_boundary is only that the 1287 1354 field_boundary will allow you to change the level of the stage height when 1288 you read in the boundary condition. This is very useful when running 1289 different tide heights in the same area as you need only to convert one 1290 boundary condition to a SWW file, ideally for tide height of 0 m 1355 you read in the boundary condition. This is very useful when running 1356 different tide heights in the same area as you need only to convert one 1357 boundary condition to a SWW file, ideally for tide height of 0 m 1291 1358 (saving disk space). Then you can use field_boundary to read this SWW file 1292 1359 and change the stage height (tide) on the fly depending on the scenario. 1293 1294 1360 """ 1295 1361 1296 1297 def __init__(self, filename, domain, 1362 ## 1363 # @brief Construct an instance of a 'field' boundary. 1364 # @param filename Name of SWW file containing stage, x and ymomentum. 1365 # @param domain Shallow water domain for which the boundary applies. 1366 # @param mean_stage Mean water level added to stage derived from SWW file. 1367 # @param time_thinning Time step thinning factor. 1368 # @param time_limit 1369 # @param boundary_polygon 1370 # @param default_boundary None or an instance of Boundary. 1371 # @param use_cache True if caching is to be used. 1372 # @param verbose True if this method is to be verbose. 1373 def __init__(self, 1374 filename, 1375 domain, 1298 1376 mean_stage=0.0, 1299 1377 time_thinning=1, 1300 1378 time_limit=None, 1301 boundary_polygon=None, 1302 default_boundary=None, 1379 boundary_polygon=None, 1380 default_boundary=None, 1303 1381 use_cache=False, 1304 1382 verbose=False): … … 1310 1388 from the sww file 1311 1389 time_thinning: Will set how many time steps from the sww file read in 1312 will be interpolated to the boundary. For example if 1390 will be interpolated to the boundary. For example if 1313 1391 the sww file has 1 second time steps and is 24 hours 1314 in length it has 86400 time steps. If you set 1315 time_thinning to 1 it will read all these steps. 1392 in length it has 86400 time steps. If you set 1393 time_thinning to 1 it will read all these steps. 1316 1394 If you set it to 100 it will read every 100th step eg 1317 1395 only 864 step. This parameter is very useful to increase 1318 the speed of a model run that you are setting up 1396 the speed of a model run that you are setting up 1319 1397 and testing. 1320 1321 default_boundary: Must be either None or an instance of a 1398 1399 default_boundary: Must be either None or an instance of a 1322 1400 class descending from class Boundary. 1323 This will be used in case model time exceeds 1401 This will be used in case model time exceeds 1324 1402 that available in the underlying data. 1325 1403 1326 1404 use_cache: 1327 1405 verbose: 1328 1329 1406 """ 1330 1407 1331 1408 # Create generic file_boundary object 1332 self.file_boundary = File_boundary(filename, domain, 1409 self.file_boundary = File_boundary(filename, 1410 domain, 1333 1411 time_thinning=time_thinning, 1334 1412 time_limit=time_limit, … … 1338 1416 verbose=verbose) 1339 1417 1340 1341 1418 # Record information from File_boundary 1342 1419 self.F = self.file_boundary.F 1343 1420 self.domain = self.file_boundary.domain 1344 1421 1345 1422 # Record mean stage 1346 1423 self.mean_stage = mean_stage 1347 1424 1348 1425 ## 1426 # @note Generate a string representation of this instance. 1349 1427 def __repr__(self): 1350 1428 return 'Field boundary' 1351 1429 1352 1430 ## 1431 # @brief Calculate 'field' boundary results. 1432 # @param vol_id 1433 # @param edge_id 1353 1434 def evaluate(self, vol_id=None, edge_id=None): 1354 1435 """Return linearly interpolated values based on domain.time … … 1356 1437 vol_id and edge_id are ignored 1357 1438 """ 1358 1439 1359 1440 # Evaluate file boundary 1360 1441 q = self.file_boundary.evaluate(vol_id, edge_id) … … 1366 1447 return q 1367 1448 1368 1369 1370 #----------------------- 1449 1450 ################################################################################ 1371 1451 # Standard forcing terms 1372 #----------------------- 1373 1452 ################################################################################ 1453 1454 ## 1455 # @brief Apply gravitational pull in the presence of bed slope. 1456 # @param domain The domain to apply gravity to. 1457 # @note Wrapper for C function gravity_c(). 1374 1458 def gravity(domain): 1375 1459 """Apply gravitational pull in the presence of bed slope … … 1377 1461 """ 1378 1462 1463 from shallow_water_ext import gravity as gravity_c 1464 1379 1465 xmom = domain.quantities['xmomentum'].explicit_update 1380 1466 ymom = domain.quantities['ymomentum'].explicit_update … … 1388 1474 x = domain.get_vertex_coordinates() 1389 1475 g = domain.g 1390 1391 1392 from shallow_water_ext import gravity as gravity_c 1393 gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6) 1394 1395 1396 1476 1477 gravity_c(g, h, z, x, xmom, ymom) #, 1.0e-6) 1478 1479 ## 1480 # @brief Apply friction to a surface (implicit). 1481 # @param domain The domain to apply Manning friction to. 1482 # @note Wrapper for C function manning_friction_c(). 1397 1483 def manning_friction_implicit(domain): 1398 """Apply (Manning) friction to water momentum 1484 """Apply (Manning) friction to water momentum 1399 1485 Wrapper for c version 1400 1486 """ 1401 1487 1402 1403 #print 'Implicit friction' 1488 from shallow_water_ext import manning_friction as manning_friction_c 1404 1489 1405 1490 xmom = domain.quantities['xmomentum'] … … 1420 1505 g = domain.g 1421 1506 1422 from shallow_water_ext import manning_friction as manning_friction_c1423 1507 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1424 1508 1425 1509 ## 1510 # @brief Apply friction to a surface (explicit). 1511 # @param domain The domain to apply Manning friction to. 1512 # @note Wrapper for C function manning_friction_c(). 1426 1513 def manning_friction_explicit(domain): 1427 """Apply (Manning) friction to water momentum 1514 """Apply (Manning) friction to water momentum 1428 1515 Wrapper for c version 1429 1516 """ 1430 1517 1431 # print 'Explicit friction'1518 from shallow_water_ext import manning_friction as manning_friction_c 1432 1519 1433 1520 xmom = domain.quantities['xmomentum'] … … 1448 1535 g = domain.g 1449 1536 1450 from shallow_water_ext import manning_friction as manning_friction_c1451 1537 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1452 1538 1453 1539 1454 1540 # FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?) 1455 # Is it still needed (30 Oct 2007)? 1541 ## 1542 # @brief Apply linear friction to a surface. 1543 # @param domain The domain to apply Manning friction to. 1544 # @note Is this still used (30 Oct 2007)? 1456 1545 def linear_friction(domain): 1457 1546 """Apply linear friction to water momentum … … 1486 1575 ymom_update[k] += S*vh[k] 1487 1576 1488 1489 1490 #--------------------------------- 1577 ################################################################################ 1491 1578 # Experimental auxiliary functions 1492 #--------------------------------- 1579 ################################################################################ 1580 1581 ## 1582 # @brief Check forcefield parameter. 1583 # @param f Object to check. 1584 # @note 'f' may be a callable object or a scalar value. 1493 1585 def check_forcefield(f): 1494 """Check that f is either 1586 """Check that force object is as expected. 1587 1588 Check that f is either: 1495 1589 1: a callable object f(t,x,y), where x and y are vectors 1496 1590 and that it returns an array or a list of same length … … 1516 1610 msg += 'not be converted into a numeric array of floats.\n' 1517 1611 msg += 'Specified function should return either list or array.' 1518 raise msg1612 raise Exception, msg 1519 1613 1520 1614 # Is this really what we want? 1521 msg = 'Return vector from function %s ' %f 1522 msg += 'must have same length as input vectors' 1523 msg += ' (type(q)=%s' % type(q) 1615 msg = 'Function %s must return vector' % str(f) 1616 assert hasattr(q, 'len'), msg 1617 1618 msg = ('Return vector from function %s must have same ' 1619 'length as input vectors' % f) 1524 1620 assert len(q) == N, msg 1525 1526 1621 else: 1527 1622 try: 1528 1623 f = float(f) 1529 1624 except: 1530 msg = 'Force field %s must be either a scalar' %f 1531 msg += ' or a vector function' 1532 raise Exception(msg) 1625 msg = ('Force field %s must be either a vector function or a ' 1626 'scalar value (coercible to float).' % str(f)) 1627 raise Exception, msg 1628 1533 1629 return f 1534 1630 1535 1631 1632 ## 1633 # Class to apply a wind stress to a domain. 1536 1634 class Wind_stress: 1537 1635 """Apply wind stress to water momentum in terms of … … 1539 1637 """ 1540 1638 1639 ## 1640 # @brief Create an instance of Wind_stress. 1641 # @param *args 1642 # @param **kwargs 1541 1643 def __init__(self, *args, **kwargs): 1542 1644 """Initialise windfield from wind speed s [m/s] … … 1591 1693 else: 1592 1694 # Assume info is in 2 keyword arguments 1593 1594 1695 if len(kwargs) == 2: 1595 1696 s = kwargs['s'] 1596 1697 phi = kwargs['phi'] 1597 1698 else: 1598 raise 'Assumes two keyword arguments: s=..., phi=....'1699 raise Exception, 'Assumes two keyword arguments: s=..., phi=....' 1599 1700 1600 1701 self.speed = check_forcefield(s) … … 1603 1704 self.const = eta_w*rho_a/rho_w 1604 1705 1605 1706 ## 1707 # @brief 'execute' this class instance. 1708 # @param domain 1606 1709 def __call__(self, domain): 1607 """Evaluate windfield based on values found in domain 1608 """ 1710 """Evaluate windfield based on values found in domain""" 1609 1711 1610 1712 from math import pi, cos, sin, sqrt … … 1613 1715 ymom_update = domain.quantities['ymomentum'].explicit_update 1614 1716 1615 N = len(domain) # number_of_triangles1717 N = len(domain) # number_of_triangles 1616 1718 t = domain.time 1617 1719 … … 1621 1723 else: 1622 1724 # Assume s is a scalar 1623 1624 1725 try: 1625 1726 s_vec = self.speed * num.ones(N, num.float) … … 1628 1729 raise msg 1629 1730 1630 1631 1731 if callable(self.phi): 1632 1732 xc = domain.get_centroid_coordinates() … … 1645 1745 1646 1746 1747 ## 1748 # @brief Assign wind field values 1749 # @param xmom_update 1750 # @param ymom_update 1751 # @param s_vec 1752 # @param phi_vec 1753 # @param const 1647 1754 def assign_windfield_values(xmom_update, ymom_update, 1648 1755 s_vec, phi_vec, const): … … 1650 1757 A c version also exists (for speed) 1651 1758 """ 1759 1652 1760 from math import pi, cos, sin, sqrt 1653 1761 … … 1670 1778 1671 1779 1672 1673 1674 1780 ## 1781 # @brief A class for a general explicit forcing term. 1675 1782 class General_forcing: 1676 1783 """General explicit forcing term for update of quantity 1677 1784 1678 1785 This is used by Inflow and Rainfall for instance 1679 1786 1680 1787 1681 1788 General_forcing(quantity_name, rate, center, radius, polygon) 1682 1789 1683 1790 domain: ANUGA computational domain 1684 quantity_name: Name of quantity to update. 1791 quantity_name: Name of quantity to update. 1685 1792 It must be a known conserved quantity. 1686 1687 rate [?/s]: Total rate of change over the specified area. 1793 1794 rate [?/s]: Total rate of change over the specified area. 1688 1795 This parameter can be either a constant or a 1689 function of time. Positive values indicate increases, 1796 function of time. Positive values indicate increases, 1690 1797 negative values indicate decreases. 1691 1798 Rate can be None at initialisation but must be specified … … 1701 1808 Either center, radius or polygon can be specified but not both. 1702 1809 If neither are specified the entire domain gets updated. 1703 All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain. 1704 1810 All coordinates to be specified in absolute UTM coordinates (x, y) assuming the zone of domain. 1811 1705 1812 Inflow or Rainfall for examples of use 1706 1813 """ … … 1709 1816 # FIXME (AnyOne) : Add various methods to allow spatial variations 1710 1817 1818 ## 1819 # @brief Create an instance of this forcing term. 1820 # @param domain 1821 # @param quantity_name 1822 # @param rate 1823 # @param center 1824 # @param radius 1825 # @param polygon 1826 # @param default_rate 1827 # @param verbose 1711 1828 def __init__(self, 1712 1829 domain, 1713 1830 quantity_name, 1714 1831 rate=0.0, 1715 center=None, radius=None, 1832 center=None, 1833 radius=None, 1716 1834 polygon=None, 1717 1835 default_rate=None, 1718 1836 verbose=False): 1719 1837 1838 from math import pi, cos, sin 1839 1720 1840 if center is None: 1721 msg = 'I got radius but no center.' 1841 msg = 'I got radius but no center.' 1722 1842 assert radius is None, msg 1723 1843 1724 1844 if radius is None: 1725 msg += 'I got center but no radius.' 1845 msg += 'I got center but no radius.' 1726 1846 assert center is None, msg 1727 1728 1729 1730 from math import pi, cos, sin1731 1847 1732 1848 self.domain = domain … … 1735 1851 self.center = ensure_numeric(center) 1736 1852 self.radius = radius 1737 self.polygon = polygon 1853 self.polygon = polygon 1738 1854 self.verbose = verbose 1739 self.value = 0.0 # Can be used to remember value at1740 # previous timestep in order to obtain rate1855 self.value = 0.0 # Can be used to remember value at 1856 # previous timestep in order to obtain rate 1741 1857 1742 1858 # Get boundary (in absolute coordinates) 1743 1859 bounding_polygon = domain.get_boundary_polygon() 1744 1860 1745 1746 1861 # Update area if applicable 1747 self.exchange_area = None 1862 self.exchange_area = None 1748 1863 if center is not None and radius is not None: 1749 1864 assert len(center) == 2 … … 1754 1869 1755 1870 # Check that circle center lies within the mesh. 1756 msg = 'Center %s specified for forcing term did not' % (str(center))1871 msg = 'Center %s specified for forcing term did not' % str(center) 1757 1872 msg += 'fall within the domain boundary.' 1758 1873 assert is_inside_polygon(center, bounding_polygon), msg … … 1762 1877 periphery_points = [] 1763 1878 for i in range(N): 1764 1765 1879 theta = 2*pi*i/100 1766 1880 1767 1881 x = center[0] + radius*cos(theta) 1768 1882 y = center[1] + radius*sin(theta) … … 1770 1884 periphery_points.append([x,y]) 1771 1885 1772 1773 1886 for point in periphery_points: 1774 msg = 'Point %s on periphery for forcing term' % (str(point))1887 msg = 'Point %s on periphery for forcing term' % str(point) 1775 1888 msg += ' did not fall within the domain boundary.' 1776 1889 assert is_inside_polygon(point, bounding_polygon), msg 1777 1890 1778 1891 if polygon is not None: 1779 1780 1892 # Check that polygon lies within the mesh. 1781 1893 for point in self.polygon: 1782 msg = 'Point %s in polygon for forcing term' % (point)1894 msg = 'Point %s in polygon for forcing term' % str(point) 1783 1895 msg += ' did not fall within the domain boundary.' 1784 1896 assert is_inside_polygon(point, bounding_polygon), msg 1785 1786 # Compute area and check that it is greater than 0 1897 1898 # Compute area and check that it is greater than 0 1787 1899 self.exchange_area = polygon_area(self.polygon) 1788 1789 msg = 'Polygon %s in forcing term' %(self.polygon) 1790 msg += ' has area = %f' %self.exchange_area 1791 assert self.exchange_area > 0.0 1792 1793 1794 1900 1901 msg = 'Polygon %s in forcing term' % str(self.polygon) 1902 msg += ' has area = %f' % self.exchange_area 1903 assert self.exchange_area > 0.0 1795 1904 1796 1905 # Pointer to update vector … … 1798 1907 1799 1908 # Determine indices in flow area 1800 N = len(domain) 1909 N = len(domain) 1801 1910 points = domain.get_centroid_coordinates(absolute=True) 1802 1911 … … 1805 1914 if self.center is not None and self.radius is not None: 1806 1915 # Inlet is circular 1807 1808 inlet_region = 'center=%s, radius=%s' %(self.center, self.radius) 1809 1916 inlet_region = 'center=%s, radius=%s' % (self.center, self.radius) 1917 1810 1918 self.exchange_indices = [] 1811 1919 for k in range(N): 1812 x, y = points[k,:] # Centroid1813 1920 x, y = points[k,:] # Centroid 1921 1814 1922 c = self.center 1815 1923 if ((x-c[0])**2+(y-c[1])**2) < self.radius**2: 1816 1924 self.exchange_indices.append(k) 1817 1818 if self.polygon is not None: 1925 1926 if self.polygon is not None: 1819 1927 # Inlet is polygon 1820 1821 inlet_region = 'polygon=%s, area=%f m^2' %(self.polygon, 1822 self.exchange_area) 1823 1928 inlet_region = 'polygon=%s, area=%f m^2' % (self.polygon, 1929 self.exchange_area) 1930 1824 1931 self.exchange_indices = inside_polygon(points, self.polygon) 1825 1826 1827 1932 1828 1933 if self.exchange_indices is not None: 1829 #print inlet_region1830 1831 1934 if len(self.exchange_indices) == 0: 1832 1935 msg = 'No triangles have been identified in ' 1833 msg += 'specified region: %s' % inlet_region1936 msg += 'specified region: %s' % inlet_region 1834 1937 raise Exception, msg 1835 1938 1836 1939 # Check and store default_rate 1837 msg = 'Keyword argument default_rate must be either None '1838 msg += 'or a function of time.\n I got %s' %(str(default_rate))1839 assert default_rate is None or \1840 type(default_rate) in [IntType, FloatType] or \1841 callable(default_rate), msg1842 1940 msg = ('Keyword argument default_rate must be either None ' 1941 'or a function of time.\nI got %s.' % str(default_rate)) 1942 assert (default_rate is None or 1943 type(default_rate) in [IntType, FloatType] or 1944 callable(default_rate)), msg 1945 1843 1946 if default_rate is not None: 1844 1845 1947 # If it is a constant, make it a function 1846 1948 if not callable(default_rate): … … 1848 1950 default_rate = lambda t: tmp 1849 1951 1850 1851 1952 # Check that default_rate is a function of one argument 1852 1953 try: … … 1856 1957 1857 1958 self.default_rate = default_rate 1858 self.default_rate_invoked = False # Flag 1859 1860 1959 self.default_rate_invoked = False # Flag 1960 1961 ## 1962 # @brief Execute this instance. 1963 # @param domain 1861 1964 def __call__(self, domain): 1862 """Apply inflow function at time specified in domain and update stage 1863 """ 1965 """Apply inflow function at time specified in domain, update stage""" 1864 1966 1865 1967 # Call virtual method allowing local modifications 1866 1867 1968 t = domain.get_time() 1868 1969 try: … … 1872 1973 except Modeltime_too_late, e: 1873 1974 if self.default_rate is None: 1874 raise Exception, e # Reraise exception1975 raise Exception, e # Reraise exception 1875 1976 else: 1876 1977 # Pass control to default rate function 1877 1978 rate = self.default_rate(t) 1878 1979 1879 1980 if self.default_rate_invoked is False: 1880 1981 # Issue warning the first time 1881 msg = '%s' %str(e)1882 msg += 'Instead I will use the default rate: %s\n'\1883 %str(self.default_rate)1884 msg += 'Note: Further warnings will be supressed'1982 msg = ('%s\n' 1983 'Instead I will use the default rate: %s\n' 1984 'Note: Further warnings will be supressed' 1985 % (str(e), str(self.default_rate))) 1885 1986 warn(msg) 1886 1987 1887 1988 # FIXME (Ole): Replace this crude flag with 1888 1989 # Python's ability to print warnings only once. 1889 1990 # See http://docs.python.org/lib/warning-filter.html 1890 1991 self.default_rate_invoked = True 1891 1892 1893 1894 1895 1992 1896 1993 if rate is None: 1897 msg = 'Attribute rate must be specified in General_forcing'1898 msg += ' or its descendants before attempting to call it'1994 msg = ('Attribute rate must be specified in General_forcing ' 1995 'or its descendants before attempting to call it') 1899 1996 raise Exception, msg 1900 1901 1997 1902 1998 # Now rate is a number 1903 1999 if self.verbose is True: 1904 print 'Rate of %s at time = %.2f = %f' %(self.quantity_name, 1905 domain.get_time(), 1906 rate) 1907 2000 print 'Rate of %s at time = %.2f = %f' % (self.quantity_name, 2001 domain.get_time(), 2002 rate) 1908 2003 1909 2004 if self.exchange_indices is None: … … 1914 2009 self.update[k] += rate 1915 2010 1916 2011 ## 2012 # @brief Update the internal rate. 2013 # @param t A callable or scalar used to set the rate. 2014 # @return The new rate. 1917 2015 def update_rate(self, t): 1918 2016 """Virtual method allowing local modifications by writing an 1919 2017 overriding version in descendant 1920 1921 """ 2018 """ 2019 1922 2020 if callable(self.rate): 1923 2021 rate = self.rate(t) … … 1927 2025 return rate 1928 2026 1929 2027 ## 2028 # @brief Get values for the specified quantity. 2029 # @param quantity_name Name of the quantity of interest. 2030 # @return The value(s) of the quantity. 2031 # @note If 'quantity_name' is None, use self.quantity_name. 1930 2032 def get_quantity_values(self, quantity_name=None): 1931 """Return values for specified quantity restricted to opening 1932 2033 """Return values for specified quantity restricted to opening 2034 1933 2035 Optionally a quantity name can be specified if values from another 1934 2036 quantity is sought 1935 2037 """ 1936 2038 1937 2039 if quantity_name is None: 1938 2040 quantity_name = self.quantity_name 1939 2041 1940 2042 q = self.domain.quantities[quantity_name] 1941 return q.get_values(location='centroids', 2043 return q.get_values(location='centroids', 1942 2044 indices=self.exchange_indices) 1943 1944 2045 2046 ## 2047 # @brief Set value for the specified quantity. 2048 # @param val The value object used to set value. 2049 # @param quantity_name Name of the quantity of interest. 2050 # @note If 'quantity_name' is None, use self.quantity_name. 1945 2051 def set_quantity_values(self, val, quantity_name=None): 1946 """Set values for specified quantity restricted to opening 1947 2052 """Set values for specified quantity restricted to opening 2053 1948 2054 Optionally a quantity name can be specified if values from another 1949 quantity is sought 2055 quantity is sought 1950 2056 """ 1951 2057 1952 2058 if quantity_name is None: 1953 2059 quantity_name = self.quantity_name 1954 1955 q = self.domain.quantities[self.quantity_name] 1956 q.set_values(val, 1957 location='centroids', 1958 indices=self.exchange_indices) 1959 1960 1961 2060 2061 q = self.domain.quantities[self.quantity_name] 2062 q.set_values(val, 2063 location='centroids', 2064 indices=self.exchange_indices) 2065 2066 2067 ## 2068 # @brief A class for rainfall forcing function. 2069 # @note Inherits from General_forcing. 1962 2070 class Rainfall(General_forcing): 1963 2071 """Class Rainfall - general 'rain over entire domain' forcing term. 1964 2072 1965 2073 Used for implementing Rainfall over the entire domain. 1966 2074 1967 2075 Current Limited to only One Gauge.. 1968 1969 Need to add Spatial Varying Capability 2076 2077 Need to add Spatial Varying Capability 1970 2078 (This module came from copying and amending the Inflow Code) 1971 2079 1972 2080 Rainfall(rain) 1973 2081 1974 domain 1975 rain [mm/s]: Total rain rate over the specified domain. 2082 domain 2083 rain [mm/s]: Total rain rate over the specified domain. 1976 2084 NOTE: Raingauge Data needs to reflect the time step. 1977 2085 IE: if Gauge is mm read at a time step, then the input … … 1980 2088 1981 2089 This parameter can be either a constant or a 1982 function of time. Positive values indicate inflow, 2090 function of time. Positive values indicate inflow, 1983 2091 negative values indicate outflow. 1984 2092 (and be used for Infiltration - Write Seperate Module) … … 1988 2096 1989 2097 polygon: Specifies a polygon to restrict the rainfall. 1990 2098 1991 2099 Examples 1992 2100 How to put them in a run File... … … 1998 2106 # input of Rainfall in mm/s 1999 2107 2000 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 2108 catchmentrainfall = Rainfall(rain=file_function('Q100_2hr_Rain.tms')) 2001 2109 # Note need path to File in String. 2002 2110 # Else assumed in same directory … … 2005 2113 """ 2006 2114 2007 2115 ## 2116 # @brief Create an instance of the class. 2117 # @param domain Domain of interest. 2118 # @param rate Total rain rate over the specified domain (mm/s). 2119 # @param center 2120 # @param radius 2121 # @param polygon Polygon to restrict rainfall. 2122 # @param default_rate 2123 # @param verbose True if this instance is to be verbose. 2008 2124 def __init__(self, 2009 2125 domain, 2010 2126 rate=0.0, 2011 center=None, radius=None, 2127 center=None, 2128 radius=None, 2012 2129 polygon=None, 2013 default_rate=None, 2130 default_rate=None, 2014 2131 verbose=False): 2015 2132 … … 2020 2137 rain = rate/1000.0 2021 2138 2022 if default_rate is not None: 2139 if default_rate is not None: 2023 2140 if callable(default_rate): 2024 2141 default_rain = lambda t: default_rate(t)/1000.0 … … 2028 2145 default_rain = None 2029 2146 2030 2031 2032 2147 # pass to parent class 2033 2148 General_forcing.__init__(self, 2034 2149 domain, 2035 2150 'stage', 2036 2151 rate=rain, 2037 center=center, radius=radius, 2152 center=center, 2153 radius=radius, 2038 2154 polygon=polygon, 2039 2155 default_rate=default_rain, 2040 2156 verbose=verbose) 2041 2157 2042 2043 2044 2045 2046 2158 2159 ## 2160 # @brief A class for inflow (rain and drain) forcing function. 2161 # @note Inherits from General_forcing. 2047 2162 class Inflow(General_forcing): 2048 2163 """Class Inflow - general 'rain and drain' forcing term. 2049 2164 2050 2165 Useful for implementing flows in and out of the domain. 2051 2166 2052 2167 Inflow(flow, center, radius, polygon) 2053 2168 2054 2169 domain 2055 rate [m^3/s]: Total flow rate over the specified area. 2170 rate [m^3/s]: Total flow rate over the specified area. 2056 2171 This parameter can be either a constant or a 2057 function of time. Positive values indicate inflow, 2172 function of time. Positive values indicate inflow, 2058 2173 negative values indicate outflow. 2059 2174 The specified flow will be divided by the area of 2060 the inflow region and then applied to update stage. 2175 the inflow region and then applied to update stage. 2061 2176 center [m]: Coordinates at center of flow point 2062 2177 radius [m]: Size of circular area … … 2064 2179 2065 2180 Either center, radius or polygon must be specified 2066 2181 2067 2182 Examples 2068 2183 2069 2184 # Constant drain at 0.003 m^3/s. 2070 2185 # The outflow area is 0.07**2*pi=0.0154 m^2 2071 # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s 2072 # 2186 # This corresponds to a rate of change of 0.003/0.0154 = 0.2 m/s 2187 # 2073 2188 Inflow((0.7, 0.4), 0.07, -0.003) 2074 2189 … … 2076 2191 # Tap turning up to a maximum inflow of 0.0142 m^3/s. 2077 2192 # The inflow area is 0.03**2*pi = 0.00283 m^2 2078 # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s 2193 # This corresponds to a rate of change of 0.0142/0.00283 = 5 m/s 2079 2194 # over the specified area 2080 2195 Inflow((0.5, 0.5), 0.03, lambda t: min(0.01*t, 0.0142)) … … 2091 2206 2092 2207 domain.forcing_terms.append(hydrograph) 2093 2094 2208 """ 2095 2209 2096 2210 ## 2211 # @brief Create an instance of the class. 2212 # @param domain Domain of interest. 2213 # @param rate Total rain rate over the specified domain (mm/s). 2214 # @param center 2215 # @param radius 2216 # @param polygon Polygon to restrict rainfall. 2217 # @param default_rate 2218 # @param verbose True if this instance is to be verbose. 2097 2219 def __init__(self, 2098 2220 domain, 2099 2221 rate=0.0, 2100 center=None, radius=None, 2222 center=None, 2223 radius=None, 2101 2224 polygon=None, 2102 2225 default_rate=None, 2103 verbose=False): 2104 2105 2226 verbose=False): 2106 2227 # Create object first to make area is available 2107 2228 General_forcing.__init__(self, … … 2109 2230 'stage', 2110 2231 rate=rate, 2111 center=center, radius=radius, 2232 center=center, 2233 radius=radius, 2112 2234 polygon=polygon, 2113 2235 default_rate=default_rate, 2114 2236 verbose=verbose) 2115 2237 2238 ## 2239 # @brief Update the instance rate. 2240 # @param t New rate object. 2116 2241 def update_rate(self, t): 2117 2242 """Virtual method allowing local modifications by writing an 2118 2243 overriding version in descendant 2119 2244 2120 This one converts m^3/s to m/s which can be added directly 2245 This one converts m^3/s to m/s which can be added directly 2121 2246 to 'stage' in ANUGA 2122 2247 """ … … 2130 2255 2131 2256 2132 2133 2134 #------------------ 2257 ################################################################################ 2135 2258 # Initialise module 2136 #------------------ 2137 2259 ################################################################################ 2138 2260 2139 2261 from anuga.utilities import compile 2140 2262 if compile.can_use_C_extension('shallow_water_ext.c'): 2141 # Underlying C implementations can be accessed 2142 2263 # Underlying C implementations can be accessed 2143 2264 from shallow_water_ext import rotate, assign_windfield_values 2144 2265 else: 2145 msg = 'C implementations could not be accessed by %s.\n ' % __file__2266 msg = 'C implementations could not be accessed by %s.\n ' % __file__ 2146 2267 msg += 'Make sure compile_all.py has been run as described in ' 2147 2268 msg += 'the ANUGA installation guide.' 2148 2269 raise Exception, msg 2149 2150 2270 2151 2271 # Optimisation with psyco … … 2160 2280 #Psyco isn't supported on 64 bit systems, but it doesn't matter 2161 2281 else: 2162 msg = 'WARNING: psyco (speedup) could not import'+\2163 ', you may want to consider installing it'2282 msg = ('WARNING: psyco (speedup) could not be imported, ' 2283 'you may want to consider installing it') 2164 2284 print msg 2165 2285 else: … … 2167 2287 psyco.bind(Domain.compute_fluxes) 2168 2288 2289 2169 2290 if __name__ == "__main__": 2170 2291 pass 2171 2172 -
branches/numpy/anuga/shallow_water/shallow_water_ext.c
r6304 r6410 879 879 } 880 880 881 // check that numpy array objects arrays are C contiguous memory 882 CHECK_C_CONTIG(h); 883 CHECK_C_CONTIG(v); 884 CHECK_C_CONTIG(x); 885 CHECK_C_CONTIG(xmom); 886 CHECK_C_CONTIG(ymom); 887 881 888 N = h -> dimensions[0]; 882 889 for (k=0; k<N; k++) { … … 937 944 } 938 945 946 // check that numpy array objects arrays are C contiguous memory 947 CHECK_C_CONTIG(w); 948 CHECK_C_CONTIG(z); 949 CHECK_C_CONTIG(uh); 950 CHECK_C_CONTIG(vh); 951 CHECK_C_CONTIG(eta); 952 CHECK_C_CONTIG(xmom); 953 CHECK_C_CONTIG(ymom); 939 954 940 955 N = w -> dimensions[0]; … … 967 982 &xmom, &ymom)) 968 983 return NULL; 984 985 // check that numpy array objects arrays are C contiguous memory 986 CHECK_C_CONTIG(w); 987 CHECK_C_CONTIG(z); 988 CHECK_C_CONTIG(uh); 989 CHECK_C_CONTIG(vh); 990 CHECK_C_CONTIG(eta); 991 CHECK_C_CONTIG(xmom); 992 CHECK_C_CONTIG(ymom); 969 993 970 994 N = w -> dimensions[0]; … … 1499 1523 } 1500 1524 1525 // check that numpy array objects arrays are C contiguous memory 1526 CHECK_C_CONTIG(surrogate_neighbours); 1527 CHECK_C_CONTIG(number_of_boundaries); 1528 CHECK_C_CONTIG(centroid_coordinates); 1529 CHECK_C_CONTIG(stage_centroid_values); 1530 CHECK_C_CONTIG(xmom_centroid_values); 1531 CHECK_C_CONTIG(ymom_centroid_values); 1532 CHECK_C_CONTIG(elevation_centroid_values); 1533 CHECK_C_CONTIG(vertex_coordinates); 1534 CHECK_C_CONTIG(stage_vertex_values); 1535 CHECK_C_CONTIG(xmom_vertex_values); 1536 CHECK_C_CONTIG(ymom_vertex_values); 1537 CHECK_C_CONTIG(elevation_vertex_values); 1538 1501 1539 // Get the safety factor beta_w, set in the config.py file. 1502 1540 // This is used in the limiting process … … 2052 2090 } 2053 2091 2054 2092 // check that numpy array objects arrays are C contiguous memory 2093 CHECK_C_CONTIG(neighbours); 2094 CHECK_C_CONTIG(neighbour_edges); 2095 CHECK_C_CONTIG(normals); 2096 CHECK_C_CONTIG(edgelengths); 2097 CHECK_C_CONTIG(radii); 2098 CHECK_C_CONTIG(areas); 2099 CHECK_C_CONTIG(tri_full_flag); 2100 CHECK_C_CONTIG(stage_edge_values); 2101 CHECK_C_CONTIG(xmom_edge_values); 2102 CHECK_C_CONTIG(ymom_edge_values); 2103 CHECK_C_CONTIG(bed_edge_values); 2104 CHECK_C_CONTIG(stage_boundary_values); 2105 CHECK_C_CONTIG(xmom_boundary_values); 2106 CHECK_C_CONTIG(ymom_boundary_values); 2107 CHECK_C_CONTIG(stage_explicit_update); 2108 CHECK_C_CONTIG(xmom_explicit_update); 2109 CHECK_C_CONTIG(ymom_explicit_update); 2110 CHECK_C_CONTIG(already_computed_flux); 2111 CHECK_C_CONTIG(max_speed_array); 2112 2055 2113 int number_of_elements = stage_edge_values -> dimensions[0]; 2056 2114 -
branches/numpy/anuga/shallow_water/test_all.py
r6304 r6410 76 76 return unittest.TestSuite(map(load, modules)) 77 77 78 ################################################################################ 79 78 80 if __name__ == '__main__': 81 79 82 suite = regressionTest() 80 83 runner = unittest.TextTestRunner() #verbosity=2) -
branches/numpy/anuga/shallow_water/test_data_manager.py
r6304 r6410 7517 7517 #print 2.0*num.transpose(ha) - stage 7518 7518 7519 ha_permutation = num.take(ha, permutation )7520 ua_permutation = num.take(ua, permutation )7521 va_permutation = num.take(va, permutation )7522 gauge_depth_permutation = num.take(gauge_depth, permutation )7519 ha_permutation = num.take(ha, permutation, axis=0) 7520 ua_permutation = num.take(ua, permutation, axis=0) 7521 va_permutation = num.take(va, permutation, axis=0) 7522 gauge_depth_permutation = num.take(gauge_depth, permutation, axis=0) 7523 7523 7524 7524 … … 7535 7535 #2.0*ha necessary because using two files with weights=1 are used 7536 7536 7537 depth_permutation = num.take(depth, permutation )7537 depth_permutation = num.take(depth, permutation, axis=0) 7538 7538 7539 7539 … … 7858 7858 # quantities written to the mux2 files subject to the permutation vector. 7859 7859 7860 ha_ref = num.take(ha0, permutation )7861 ua_ref = num.take(ua0, permutation )7862 va_ref = num.take(va0, permutation )7863 7864 gauge_depth_ref = num.take(gauge_depth, permutation )7860 ha_ref = num.take(ha0, permutation, axis=0) 7861 ua_ref = num.take(ua0, permutation, axis=0) 7862 va_ref = num.take(va0, permutation, axis=0) 7863 7864 gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) 7865 7865 7866 7866 assert num.allclose(num.transpose(ha_ref)+tide, stage0) # Meters … … 7966 7966 # quantities written to the mux2 files subject to the permutation vector. 7967 7967 7968 ha_ref = num.take(ha1, permutation )7969 ua_ref = num.take(ua1, permutation )7970 va_ref = num.take(va1, permutation )7971 7972 gauge_depth_ref = num.take(gauge_depth, permutation )7968 ha_ref = num.take(ha1, permutation, axis=0) 7969 ua_ref = num.take(ua1, permutation, axis=0) 7970 va_ref = num.take(va1, permutation, axis=0) 7971 7972 gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) 7973 7973 7974 7974 … … 8074 8074 # quantities written to the mux2 files subject to the permutation vector. 8075 8075 8076 ha_ref = weights[0]*num.take(ha0, permutation) + weights[1]*num.take(ha1, permutation) 8077 ua_ref = weights[0]*num.take(ua0, permutation) + weights[1]*num.take(ua1, permutation) 8078 va_ref = weights[0]*num.take(va0, permutation) + weights[1]*num.take(va1, permutation) 8079 8080 gauge_depth_ref = num.take(gauge_depth, permutation) 8076 ha_ref = (weights[0]*num.take(ha0, permutation, axis=0) 8077 + weights[1]*num.take(ha1, permutation, axis=0)) 8078 ua_ref = (weights[0]*num.take(ua0, permutation, axis=0) 8079 + weights[1]*num.take(ua1, permutation, axis=0)) 8080 va_ref = (weights[0]*num.take(va0, permutation, axis=0) 8081 + weights[1]*num.take(va1, permutation, axis=0)) 8082 8083 gauge_depth_ref = num.take(gauge_depth, permutation, axis=0) 8081 8084 8082 8085 … … 10588 10591 h = stage-z 10589 10592 for i in range(len(stage)): 10590 if h[i] == 0.0:10591 assert xmomentum[i] == 0.010592 assert ymomentum[i] == 0.010593 if num.alltrue(h[i] == 0.0): 10594 assert num.alltrue(xmomentum[i] == 0.0) 10595 assert num.alltrue(ymomentum[i] == 0.0) 10593 10596 else: 10594 assert h[i] >= domain.minimum_storable_height 10597 msg = ('h[i]=\n%s\ndomain.minimum_storable_height=%s' 10598 % (str(h[i]), str(domain.minimum_storable_height))) 10599 assert num.alltrue(h[i] >= domain.minimum_storable_height), msg 10595 10600 10596 10601 fid.close() … … 11329 11334 11330 11335 #------------------------------------------------------------- 11336 11331 11337 if __name__ == "__main__": 11332 11333 11338 suite = unittest.makeSuite(Test_Data_Manager,'test') 11334 #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_sts') 11335 #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo') 11336 #suite = unittest.makeSuite(Test_Data_Manager,'covered_') 11337 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources') 11338 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources') 11339 #suite = unittest.makeSuite(Test_Data_Manager,'test_get_maximum_inundation') 11339 11340 11340 11341 # FIXME (Ole): This is the test that fails under Windows 11341 11342 #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2') 11342 11343 #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV') 11343 11344 11344 11345 11345 if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V': -
branches/numpy/anuga/shallow_water/test_eq.py
r6304 r6410 60 60 61 61 #------------------------------------------------------------- 62 62 63 if __name__ == "__main__": 63 64 suite = unittest.makeSuite(Test_eq,'test_Okada_func') -
branches/numpy/anuga/shallow_water/test_shallow_water_domain.py
r6304 r6410 199 199 200 200 return 17.7 201 202 203 def scalar_func_list(t,x,y): 204 """Function that returns a scalar. 205 Used to test error message when numeric array is expected 206 """ 207 208 return [17.7] 201 209 202 210 … … 2462 2470 2463 2471 try: 2464 domain.forcing_terms.append(Wind_stress(s = scalar_func ,2472 domain.forcing_terms.append(Wind_stress(s = scalar_func_list, 2465 2473 phi = angle)) 2466 2474 except AssertionError: … … 2474 2482 domain.forcing_terms.append(Wind_stress(s = speed, 2475 2483 phi = scalar_func)) 2476 except AssertionError:2484 except Exception: 2477 2485 pass 2478 2486 else: … … 5582 5590 #print take(cv2, (0,3,8)) 5583 5591 5584 assert num.allclose( num.take(cv1, (0,8,16)), num.take(cv2, (0,3,8))) #Diag 5585 assert num.allclose( num.take(cv1, (0,6,12)), num.take(cv2, (0,1,4))) #Bottom 5586 assert num.allclose( num.take(cv1, (12,14,16)), num.take(cv2, (4,6,8))) #RHS 5592 assert num.allclose(num.take(cv1, (0,8,16), axis=0), 5593 num.take(cv2, (0,3,8), axis=0)) #Diag 5594 assert num.allclose(num.take(cv1, (0,6,12), axis=0), 5595 num.take(cv2, (0,1,4), axis=0)) #Bottom 5596 assert num.allclose(num.take(cv1, (12,14,16), axis=0), 5597 num.take(cv2, (4,6,8), axis=0)) #RHS 5587 5598 5588 5599 #Cleanup … … 5690 5701 5691 5702 #print points[0], points[5], points[10], points[15] 5692 msg = ('value was %s,\nshould be [[0,0], [1.0/3, 1.0/3], ' 5693 '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15]))) 5694 assert num.allclose(num.take(points, [0,5,10,15]), 5703 msg = ('value was\n%s\nshould be\n' 5704 '[[0,0], [1.0/3, 1.0/3],\n' 5705 '[2.0/3, 2.0/3], [1,1]]' 5706 % str(num.take(points, [0,5,10,15], axis=0))) 5707 assert num.allclose(num.take(points, [0,5,10,15], axis=0), 5695 5708 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg 5696 5709 … … 5873 5886 5874 5887 #print points[0], points[5], points[10], points[15] 5875 msg = ('values was %s,\nshould be [[0,0], [1.0/3, 1.0/3], ' 5876 '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15]))) 5877 assert num.allclose(num.take(points, [0,5,10,15]), 5888 msg = ('values was\n%s\nshould be\n' 5889 '[[0,0], [1.0/3, 1.0/3],\n' 5890 '[2.0/3, 2.0/3], [1,1]]' 5891 % str(num.take(points, [0,5,10,15], axis=0))) 5892 assert num.allclose(num.take(points, [0,5,10,15], axis=0), 5878 5893 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg 5879 5894 … … 6048 6063 6049 6064 #print points[0], points[5], points[10], points[15] 6050 msg = ('value was %s,\nshould be [[0,0], [1.0/3, 1.0/3], ' 6051 '[2.0/3, 2.0/3], [1,1]]' % str(num.take(points, [0,5,10,15]))) 6052 assert num.allclose(num.take(points, [0,5,10,15]), 6065 msg = ('value was\n%s\nshould be\n' 6066 '[[0,0], [1.0/3, 1.0/3],\n' 6067 '[2.0/3, 2.0/3], [1,1]]' 6068 % str(num.take(points, [0,5,10,15], axis=0))) 6069 assert num.allclose(num.take(points, [0,5,10,15], axis=0), 6053 6070 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]]), msg 6054 6071 … … 6626 6643 6627 6644 if __name__ == "__main__": 6628 6629 suite = unittest.makeSuite(Test_Shallow_Water, 'test') 6630 #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve') 6631 #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g') 6632 #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain') 6633 #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters') 6634 #suite = unittest.makeSuite(Test_Shallow_Water,'test_inflow_outflow_conservation') 6635 #suite = unittest.makeSuite(Test_Shallow_Water,'test_outflow_conservation_problem_temp') 6636 6637 6638 6645 #suite = unittest.makeSuite(Test_Shallow_Water, 'test') 6646 suite = unittest.makeSuite(Test_Shallow_Water, 'test_get_maximum_inundation_from_sww') 6639 6647 runner = unittest.TextTestRunner(verbosity=1) 6640 6648 runner.run(suite) -
branches/numpy/anuga/shallow_water/test_smf.py
r6304 r6410 132 132 133 133 #------------------------------------------------------------- 134 134 135 if __name__ == "__main__": 135 #suite = unittest.makeSuite(Test_smf,'test_Double_gaussian')136 136 suite = unittest.makeSuite(Test_smf,'test') 137 137 runner = unittest.TextTestRunner() -
branches/numpy/anuga/shallow_water/test_system.py
r6304 r6410 187 187 188 188 #------------------------------------------------------------- 189 189 190 if __name__ == "__main__": 190 191 suite = unittest.makeSuite(Test_system,'test') 191 #suite = unittest.makeSuite(Test_system,'test_boundary_timeII')192 192 runner = unittest.TextTestRunner() 193 193 runner.run(suite) -
branches/numpy/anuga/shallow_water/test_tsunami_okada.py
r6304 r6410 291 291 292 292 #------------------------------------------------------------- 293 293 294 if __name__ == "__main__": 294 295 suite = unittest.makeSuite(Test_eq,'test') -
branches/numpy/anuga/shallow_water/tsunami_okada.py
r6304 r6410 215 215 yrec=x 216 216 for i in range(0,len(zrec[0])): 217 if zrec[0][i]==yrec and zrec[1][i]==xrec: 217 if (num.alltrue(zrec[0][i]==yrec) and 218 num.alltrue(zrec[1][i]==xrec)): 218 219 Z=zrec[2][i] 219 220 Z=0.001*Z -
branches/numpy/anuga/utilities/polygon.py
r6360 r6410 1 1 #!/usr/bin/env python 2 """Polygon manipulations 3 """ 2 3 """Polygon manipulations""" 4 4 5 5 import numpy as num … … 11 11 12 12 13 ## 14 # @brief Determine whether a point is on a line segment. 15 # @param point (x, y) of point in question (tuple, list or array). 16 # @param line ((x1,y1), (x2,y2)) for line (tuple, list or array). 17 # @param rtol Relative error for 'close'. 18 # @param atol Absolute error for 'close'. 19 # @return True or False. 13 20 def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8): 14 21 """Determine whether a point is on a line segment 15 22 16 Input: 23 Input: 17 24 point is given by [x, y] 18 25 line is given by [x0, y0], [x1, y1]] or … … 21 28 Output: 22 29 23 Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding. 30 Note: Line can be degenerate and function still works to discern coinciding 31 points from non-coinciding. 24 32 """ 25 33 … … 31 39 line[1,0], line[1,1], 32 40 rtol, atol) 33 41 34 42 return bool(res) 35 43 … … 42 50 # result functions for possible states 43 51 def lines_dont_coincide(p0,p1,p2,p3): return (3, None) 44 def lines_0_fully_included_in_1(p0,p1,p2,p3): return (2, num.array([p0,p1])) 45 def lines_1_fully_included_in_0(p0,p1,p2,p3): return (2, num.array([p2,p3])) 46 def lines_overlap_same_direction(p0,p1,p2,p3): return (2, num.array([p0,p3])) 47 def lines_overlap_same_direction2(p0,p1,p2,p3): return (2, num.array([p2,p1])) 48 def lines_overlap_opposite_direction(p0,p1,p2,p3): return (2, num.array([p0,p2])) 49 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, num.array([p3,p1])) 52 def lines_0_fully_included_in_1(p0,p1,p2,p3): return (2, 53 num.array([p0,p1])) 54 def lines_1_fully_included_in_0(p0,p1,p2,p3): return (2, 55 num.array([p2,p3])) 56 def lines_overlap_same_direction(p0,p1,p2,p3): return (2, 57 num.array([p0,p3])) 58 def lines_overlap_same_direction2(p0,p1,p2,p3): return (2, 59 num.array([p2,p1])) 60 def lines_overlap_opposite_direction(p0,p1,p2,p3): return (2, 61 num.array([p0,p2])) 62 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, 63 num.array([p3,p1])) 50 64 51 65 # this function called when an impossible state is found 52 def lines_error(p1, p2, p3, p4): raise RuntimeError, ("INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s" % 53 (str(p1), str(p2), str(p3), str(p4))) 66 def lines_error(p1, p2, p3, p4): 67 raise RuntimeError, ('INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s' 68 % (str(p1), str(p2), str(p3), str(p4))) 54 69 55 70 # 0s1 0e1 1s0 1e0 # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0 … … 72 87 } 73 88 89 ## 90 # @brief Finds intersection point of two line segments. 91 # @param line0 First line ((x1,y1), (x2,y2)). 92 # @param line1 Second line ((x1,y1), (x2,y2)). 93 # @param rtol Relative error for 'close'. 94 # @param atol Absolute error for 'close'. 95 # @return (status, value) where: 96 # status = 0 - no intersection, value set to None 97 # 1 - intersection found, value=(x,y) 98 # 2 - lines collienar, overlap, value=overlap segment 99 # 3 - lines collinear, no overlap, value is None 100 # 4 - lines parallel, value is None 74 101 def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8): 75 """Returns intersecting point between two line segments or None 76 (if parallel or no intersection is found). 77 78 However, if parallel lines coincide partly (i.e. shara a common segment, 102 """Returns intersecting point between two line segments. 103 104 However, if parallel lines coincide partly (i.e. share a common segment), 79 105 the line segment where lines coincide is returned 80 81 106 82 107 Inputs: … … 85 110 corresponding to a point. 86 111 87 88 112 Output: 89 status, value 90 91 where status and value is interpreted as follows 92 113 status, value - where status and value is interpreted as follows: 93 114 status == 0: no intersection, value set to None. 94 115 status == 1: intersection point found and returned in value as [x,y]. 95 status == 2: Collinear overlapping lines found. Value takes the form [[x0,y0], [x1,y1]]. 116 status == 2: Collinear overlapping lines found. 117 Value takes the form [[x0,y0], [x1,y1]]. 96 118 status == 3: Collinear non-overlapping lines. Value set to None. 97 status == 4: Lines are parallel with a fixed distance apart. Value set to None. 98 119 status == 4: Lines are parallel. Value set to None. 99 120 """ 100 121 … … 102 123 103 124 line0 = ensure_numeric(line0, num.float) 104 line1 = ensure_numeric(line1, num.float) 125 line1 = ensure_numeric(line1, num.float) 105 126 106 127 x0 = line0[0,0]; y0 = line0[0,1] … … 113 134 u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2) 114 135 u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) 115 136 116 137 if num.allclose(denom, 0.0, rtol=rtol, atol=atol): 117 138 # Lines are parallel - check if they are collinear … … 123 144 point_on_line([x3, y3], line0, rtol=rtol, atol=atol)) 124 145 125 return collinear_result[state_tuple]([x0,y0],[x1,y1],[x2,y2],[x3,y3]) 146 return collinear_result[state_tuple]([x0,y0], [x1,y1], 147 [x2,y2], [x3,y3]) 126 148 else: 127 149 # Lines are parallel but aren't collinear 128 return 4, None #FIXME (Ole): Add distance here instead of None 150 return 4, None #FIXME (Ole): Add distance here instead of None 129 151 else: 130 152 # Lines are not parallel, check if they intersect 131 153 u0 = u0/denom 132 u1 = u1/denom 154 u1 = u1/denom 133 155 134 156 x = x0 + u0*(x1-x0) … … 137 159 # Sanity check - can be removed to speed up if needed 138 160 assert num.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol) 139 assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol) 161 assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol) 140 162 141 163 # Check if point found lies within given line segments 142 if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 164 if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 143 165 # We have intersection 144 166 return 1, num.array([x, y]) … … 147 169 return 0, None 148 170 149 171 ## 172 # @brief Finds intersection point of two line segments. 173 # @param line0 First line ((x1,y1), (x2,y2)). 174 # @param line1 Second line ((x1,y1), (x2,y2)). 175 # @return (status, value) where: 176 # status = 0 - no intersection, value set to None 177 # 1 - intersection found, value=(x,y) 178 # 2 - lines collienar, overlap, value=overlap segment 179 # 3 - lines collinear, no overlap, value is None 180 # 4 - lines parallel, value is None 181 # @note Wrapper for C function. 150 182 def NEW_C_intersection(line0, line1): 151 #FIXME(Ole): To write in C 152 """Returns intersecting point between two line segments or None 153 (if parallel or no intersection is found). 154 155 However, if parallel lines coincide partly (i.e. shara a common segment, 183 """Returns intersecting point between two line segments. 184 185 However, if parallel lines coincide partly (i.e. share a common segment), 156 186 the line segment where lines coincide is returned 157 158 187 159 188 Inputs: 160 189 line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]] 161 A line can also be a 2x2 num ericarray with each row190 A line can also be a 2x2 numpy array with each row 162 191 corresponding to a point. 163 192 164 165 193 Output: 166 status, value 167 168 where status is interpreted as follows 169 170 status == 0: no intersection with value set to None 171 status == 1: One intersection point found and returned in value as [x,y] 172 status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]] 173 status == 3: Lines would coincide but only if extended. Value set to None 174 status == 4: Lines are parallel with a fixed distance apart. Value set to None. 175 176 """ 177 194 status, value - where status and value is interpreted as follows: 195 status == 0: no intersection, value set to None. 196 status == 1: intersection point found and returned in value as [x,y]. 197 status == 2: Collinear overlapping lines found. 198 Value takes the form [[x0,y0], [x1,y1]]. 199 status == 3: Collinear non-overlapping lines. Value set to None. 200 status == 4: Lines are parallel. Value set to None. 201 """ 178 202 179 203 line0 = ensure_numeric(line0, num.float) 180 line1 = ensure_numeric(line1, num.float) 204 line1 = ensure_numeric(line1, num.float) 181 205 182 206 status, value = _intersection(line0[0,0], line0[0,1], … … 187 211 return status, value 188 212 189 190 191 213 ## 214 # @brief Determine if one point is inside a polygon. 215 # @param point The point of interest. 216 # @param polygon The polygon to test inclusion in. 217 # @param closed True if points on boundary are considered 'inside' polygon. 218 # @param verbose True if this function is to be verbose. 219 # @return True if point is inside the polygon. 220 # @note Uses inside_polygon() to do the work. 221 # @note Raises Exception if more than one point supplied. 192 222 def is_inside_polygon(point, polygon, closed=True, verbose=False): 193 223 """Determine if one point is inside a polygon … … 195 225 See inside_polygon for more details 196 226 """ 197 198 # print 'polygon.py: 198, is_inside_polygon: point=%s' % str(point)199 # print 'polygon.py: 198, is_inside_polygon: polygon=%s' % str(polygon)200 227 201 228 indices = inside_polygon(point, polygon, closed, verbose) … … 208 235 msg = 'is_inside_polygon must be invoked with one point only' 209 236 raise Exception, msg 210 211 237 238 ## 239 # @brief Determine which of a set of points are inside a polygon. 240 # @param points A set of points (tuple, list or array). 241 # @param polygon A set of points defining a polygon (tuple, list or array). 242 # @param closed True if points on boundary are considered 'inside' polygon. 243 # @param verbose True if this function is to be verbose. 244 # @return A list of indices of points inside the polygon. 212 245 def inside_polygon(points, polygon, closed=True, verbose=False): 213 246 """Determine points inside a polygon 214 247 215 248 Functions inside_polygon and outside_polygon have been defined in 216 terms af separate_by_polygon which will put all inside indices in249 terms of separate_by_polygon which will put all inside indices in 217 250 the first part of the indices array and outside indices in the last 218 251 … … 222 255 a list or a numeric array 223 256 """ 224 225 #if verbose: print 'Checking input to inside_polygon'226 257 227 258 try: … … 232 263 # If this fails it is going to be because the points can't be 233 264 # converted to a numeric array. 234 msg = 'Points could not be converted to numeric array' 265 msg = 'Points could not be converted to numeric array' 235 266 raise Exception, msg 236 267 … … 242 273 # If this fails it is going to be because the points can't be 243 274 # converted to a numeric array. 244 msg = 'Polygon %s could not be converted to numeric array' %(str(polygon)) 275 msg = ('Polygon %s could not be converted to numeric array' 276 % (str(polygon))) 245 277 raise Exception, msg 246 278 … … 253 285 verbose=verbose) 254 286 255 # print 'polygon.py: 255, inside_polygon: indices=%s' % str(indices)256 # print 'polygon.py: 255, inside_polygon: count=%s' % str(count)257 287 # Return indices of points inside polygon 258 288 return indices[:count] 259 289 260 261 290 ## 291 # @brief Determine if one point is outside a polygon. 292 # @param point The point of interest. 293 # @param polygon The polygon to test inclusion in. 294 # @param closed True if points on boundary are considered 'inside' polygon. 295 # @param verbose True if this function is to be verbose. 296 # @return True if point is outside the polygon. 297 # @note Uses inside_polygon() to do the work. 262 298 def is_outside_polygon(point, polygon, closed=True, verbose=False, 263 299 points_geo_ref=None, polygon_geo_ref=None): … … 268 304 269 305 indices = outside_polygon(point, polygon, closed, verbose) 270 #points_geo_ref, polygon_geo_ref)271 306 272 307 if indices.shape[0] == 1: … … 277 312 msg = 'is_outside_polygon must be invoked with one point only' 278 313 raise Exception, msg 279 280 314 315 ## 316 # @brief Determine which of a set of points are outside a polygon. 317 # @param points A set of points (tuple, list or array). 318 # @param polygon A set of points defining a polygon (tuple, list or array). 319 # @param closed True if points on boundary are considered 'inside' polygon. 320 # @param verbose True if this function is to be verbose. 321 # @return A list of indices of points outside the polygon. 281 322 def outside_polygon(points, polygon, closed = True, verbose = False): 282 323 """Determine points outside a polygon 283 324 284 325 Functions inside_polygon and outside_polygon have been defined in 285 terms af separate_by_polygon which will put all inside indices in326 terms of separate_by_polygon which will put all inside indices in 286 327 the first part of the indices array and outside indices in the last 287 328 … … 289 330 """ 290 331 291 #if verbose: print 'Checking input to outside_polygon'292 332 try: 293 333 points = ensure_numeric(points, num.float) … … 306 346 raise Exception, msg 307 347 308 309 348 if len(points.shape) == 1: 310 349 # Only one point was passed in. Convert to array of points … … 321 360 else: 322 361 return indices[count:][::-1] #return reversed 323 324 325 def in_and_outside_polygon(points, polygon, closed = True, verbose = False): 362 363 ## 364 # @brief Separate a list of points into two sets inside+outside a polygon. 365 # @param points A set of points (tuple, list or array). 366 # @param polygon A set of points defining a polygon (tuple, list or array). 367 # @param closed True if points on boundary are considered 'inside' polygon. 368 # @param verbose True if this function is to be verbose. 369 # @return A tuple (in, out) of point indices for poinst inside amd outside. 370 def in_and_outside_polygon(points, polygon, closed=True, verbose=False): 326 371 """Determine points inside and outside a polygon 327 372 328 373