Changeset 5729
- Timestamp:
- Sep 4, 2008, 2:45:27 PM (16 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r5666 r5729 369 369 370 370 def get_quantity(self, name, location='vertices', indices = None): 371 """Get quantity object.371 """Get pointer to quantity object. 372 372 373 373 name: Name of quantity -
anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r5321 r5729 1113 1113 1114 1114 1115 1115 1116 1116 1117 class Triangle_intersection: … … 1147 1148 1148 1149 return s 1150 1151 def segment_midpoints(segments): 1152 """Calculate midpoints of all segments 1153 1154 Inputs: 1155 segments: List of instances of class Segment 1156 1157 Ouputs: 1158 midpoints: List of points 1159 """ 1160 1161 midpoints = [] 1162 msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection' 1163 for segment in segments: 1164 assert isinstance(segment, Triangle_intersection), msg 1165 1166 midpoint = sum(array(segment.segment))/2 1167 midpoints.append(midpoint) 1168 1169 return midpoints 1170 1171 1172 -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
r5522 r5729 753 753 verbose=False, 754 754 use_cache=False): 755 # FIXME: Use this function for the time being. Later move code in here 755 """ Set values based on geo referenced geospatial data object. 756 """ 756 757 757 758 points = geospatial_data.get_data_points(absolute=False) 758 759 values = geospatial_data.get_attributes() 759 760 data_georef = geospatial_data.get_geo_reference() 760 761 761 762 762 … … 1007 1007 1008 1008 def get_interpolated_values(self, interpolation_points): 1009 1010 # Interpolation object based on internal (discontinuous triangles) 1009 """ Get values at interpolation points 1010 1011 If interpolation points have been given previously, the 1012 associated matrices will be reused to save time. 1013 1014 The argument interpolation points must be given as either a list of absolute UTM coordinates or 1015 a geospatial data object. 1016 """ 1017 1018 # FIXME (Ole): Could do with an input check (should be generalised and used widely) 1019 # That will check that interpolation points is either a list of points, Nx2 array, or geospatial 1020 1021 # Ensure points are converted to coordinates relative to mesh origin 1022 # FIXME (Ole): This could all be refactored using the 'change_points_geo_ref' method 1023 # of Class geo_reference. The purpose is to make interpolation points relative 1024 # to the mesh origin. 1025 # 1026 # Speed is also a consideration here. 1027 1028 if isinstance(interpolation_points, Geospatial_data): 1029 # Ensure interpolation points are in absolute UTM coordinates 1030 interpolation_points = interpolation_points.get_data_points(absolute=True) 1031 1032 # Reconcile interpolation points with georeference of domain 1033 interpolation_points = self.domain.geo_reference.get_relative(interpolation_points) 1034 interpolation_points = ensure_numeric(interpolation_points) 1035 1036 # Get internal (discontinuous) triangles for use with the interpolation object. 1011 1037 x, y, vertex_values, triangles = self.get_vertex_values(xy=True, 1012 1038 smooth=False) 1013 # FIXME : This concat should roll into get_vertex_values1039 # FIXME (Ole): This concat should roll into get_vertex_values 1014 1040 vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]), 1015 1041 axis=1) … … 1020 1046 I = self.interpolation_object 1021 1047 1022 if allclose(interpolation_points, I._point_coordinates): 1023 can_reuse = True 1048 if allclose(interpolation_points.shape, 1049 I._point_coordinates.shape): 1050 if allclose(interpolation_points, I._point_coordinates): 1051 can_reuse = True 1024 1052 1025 1053 … … 1047 1075 1048 1076 return X, Compatible list, Numeric array (see below) 1049 interpolation_points: List of x, y coordinates where value is 1050 sought (using interpolation). If points are given, values of 1051 location and indices are ignored 1052 1053 location: Where values are to be stored. 1054 Permissible options are: vertices, edges, centroids 1055 and unique vertices. Default is 'vertices' 1077 1078 Inputs: 1079 interpolation_points: List of x, y coordinates where value is 1080 sought (using interpolation). If points 1081 are given, values of location and indices 1082 are ignored. Assume either absolute UTM 1083 coordinates or geospatial data object. 1084 1085 location: Where values are to be stored. 1086 Permissible options are: vertices, edges, centroids 1087 and unique vertices. Default is 'vertices' 1056 1088 1057 1089 … … 1075 1107 internal ordering. 1076 1108 """ 1109 1077 1110 from Numeric import take 1078 1111 … … 1086 1119 1087 1120 # FIXME (Ole): Consider deprecating 'edges' - but not if it is used 1088 # elsewhere in ANUGA. 1121 # elsewhere in ANUGA. 1122 # Edges have already been deprecated in set_values, see changeset:5521, 1123 # but *might* be useful in get_values. Any thoughts anyone? 1124 1089 1125 if location not in ['vertices', 'centroids', 'edges', 1090 1126 'unique vertices']: -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r5521 r5729 2354 2354 2355 2355 2356 def test_get_interpolated_values_with_georef(self): 2357 2358 zone = 56 2359 xllcorner = 308500 2360 yllcorner = 6189000 2361 a = [0.0, 0.0] 2362 b = [0.0, 2.0] 2363 c = [2.0,0.0] 2364 d = [0.0, 4.0] 2365 e = [2.0, 2.0] 2366 f = [4.0,0.0] 2367 2368 points = [a, b, c, d, e, f] 2369 #bac, bce, ecf, dbe 2370 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2371 2372 domain = Domain(points, vertices, 2373 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2374 2375 quantity = Quantity(domain) 2376 quantity.set_values(lambda x, y: x+2*y) #2 4 4 6 2377 2378 #First pick one point (and turn it into absolute coordinates) 2379 x, y = 2.0/3, 8.0/3 2380 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2381 assert allclose(v, 6) 2382 2383 2384 # Then another to test that algorithm won't blindly 2385 # reuse interpolation matrix 2386 x, y = 4.0/3, 4.0/3 2387 v = quantity.get_values(interpolation_points = [[x+xllcorner,y+yllcorner]]) 2388 assert allclose(v, 4) 2389 2390 # Try two points 2391 pts = [[2.0/3 + xllcorner, 8.0/3 + yllcorner], 2392 [4.0/3 + xllcorner, 4.0/3 + yllcorner]] 2393 v = quantity.get_values(interpolation_points=pts) 2394 assert allclose(v, [6, 4]) 2395 2396 # Test it using the geospatial data format with absolute input points and default georef 2397 pts = Geospatial_data(data_points=pts) 2398 v = quantity.get_values(interpolation_points=pts) 2399 assert allclose(v, [6, 4]) 2400 2401 2402 # Test it using the geospatial data format with relative input points 2403 pts = Geospatial_data(data_points=[[2.0/3, 8.0/3], [4.0/3, 4.0/3]], 2404 geo_reference=Geo_reference(zone,xllcorner,yllcorner)) 2405 v = quantity.get_values(interpolation_points=pts) 2406 assert allclose(v, [6, 4]) 2407 2408 2409 2356 2410 2357 2411 def test_getting_some_vertex_values(self): -
anuga_core/source/anuga/geospatial_data/geospatial_data.py
r5421 r5729 1543 1543 alpha = alpha) 1544 1544 1545 points_geo=domain.geo_reference.change_points_geo_ref(points) 1545 1546 # Convert points to geospatial data for use with get_values below 1547 points_geo = Geospatial_data(points, domain.geo_reference) 1548 1546 1549 #returns the predicted elevation of the points that were "split" out 1547 1550 #of the original data set for one particular alpha -
anuga_core/source/anuga/geospatial_data/test_geospatial_data.py
r5203 r5729 1892 1892 os.remove(filename) 1893 1893 1894 #print value, alpha1894 # print value, alpha 1895 1895 assert (alpha==0.01) 1896 1896 … … 1976 1976 1977 1977 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long') 1978 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')1979 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')1978 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter') 1979 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1') 1980 1980 suite = unittest.makeSuite(Test_Geospatial_data, 'test') 1981 1981 runner = unittest.TextTestRunner() #verbosity=2) -
anuga_core/source/anuga/shallow_water/data_manager.py
r5726 r5729 100 100 from anuga.abstract_2d_finite_volumes.util import get_revision_number, \ 101 101 remove_lone_verts, sww2timeseries, get_centroid_values 102 103 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints 102 104 from anuga.load_mesh.loadASCII import export_mesh_file 103 105 from anuga.utilities.polygon import intersection … … 6370 6372 segments = mesh.get_intersecting_segments(polyline, verbose=verbose) 6371 6373 6372 6373 # Then store for each triangle the length of the intersecting segment(s), 6374 # right hand normal(s) and midpoints as interpolation_points. 6375 interpolation_points = [] 6376 for segment in segments: 6377 midpoint = sum(array(segment.segment))/2 6378 interpolation_points.append(midpoint) 6379 6380 6374 # Get midpoints 6375 interpolation_points = segment_midpoints(segments) 6376 6377 # Interpolate 6381 6378 if verbose: 6382 6379 print 'Interpolating - ', 6383 6380 print 'total number of interpolation points = %d'\ 6384 6381 %len(interpolation_points) 6385 6386 6387 6382 6388 6383 I = Interpolation_function(time, … … 6445 6440 for t in time: 6446 6441 total_flow=0 6447 for i , p in enumerate(interpolation_points):6442 for i in range(len(interpolation_points)): 6448 6443 elevation, stage, uh, vh = interpolation_function(t, point_id=i) 6449 6444 normal = segments[i].normal -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r5672 r5729 83 83 from Numeric import compress, arange 84 84 85 85 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints 86 86 from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain 87 87 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ … … 374 374 return self.get_quantity('elevation').\ 375 375 get_maximum_location(indices=wet_elements) 376 377 378 379 380 def get_flow_through_cross_section(self, polyline, 381 verbose=False): 382 """Get the total flow through an arbitrary poly line. 383 384 This is a run-time equivalent of the function with same name in data_manager.py 385 386 Input: 387 polyline: Representation of desired cross section - it may contain 388 multiple sections allowing for complex shapes. Assume 389 absolute UTM coordinates. 390 Format [[x0, y0], [x1, y1], ...] 391 392 Output: 393 Q: Total flow [m^3/s] across given segments. 394 395 396 """ 397 398 # Adjust polyline to mesh spatial origin 399 polyline = self.geo_reference.get_relative(polyline) 400 401 # Find all intersections and associated triangles. 402 segments = self.get_intersecting_segments(polyline, verbose=verbose) 403 404 msg = 'No segments found' 405 assert len(segments) > 0, msg 406 407 # Get midpoints 408 midpoints = segment_midpoints(segments) 409 410 # FIXME (Ole): HACK - need to make midpoints Geospatial instances 411 midpoints = self.geo_reference.get_absolute(midpoints) 412 413 # Compute flow 414 if verbose: print 'Computing flow through specified cross section' 415 416 # Get interpolated values 417 xmomentum = self.get_quantity('xmomentum') 418 ymomentum = self.get_quantity('ymomentum') 419 420 uh = xmomentum.get_values(interpolation_points=midpoints) 421 vh = ymomentum.get_values(interpolation_points=midpoints) 422 423 # Compute and sum flows across each segment 424 total_flow=0 425 for i in range(len(uh)): 426 427 # Inner product of momentum vector with segment normal [m^2/s] 428 normal = segments[i].normal 429 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1] 430 431 # Flow across this segment [m^3/s] 432 segment_flow = normal_momentum*segments[i].length 433 434 # Accumulate 435 total_flow += segment_flow 436 437 438 return total_flow 439 376 440 377 441 def check_integrity(self): -
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5726 r5729 10073 10073 u = 2 m/s 10074 10074 h = 1 m 10075 w = 5m (width of channel)10076 10077 q = u*h*w = 10m^3/s10075 w = 3 m (width of channel) 10076 10077 q = u*h*w = 6 m^3/s 10078 10078 10079 10079 #---------- First run without geo referencing … … 10112 10112 10113 10113 Br = Reflective_boundary(domain) # Side walls 10114 Bd = Dirichlet_boundary([h, uh, 0]) # 2 m/s across the 5m inlet:10114 Bd = Dirichlet_boundary([h, uh, 0]) # 2 m/s across the 3 m inlet: 10115 10115 10116 10116 … … 10199 10199 u = 2 m/s 10200 10200 h = 2 m 10201 w = 5m (width of channel)10202 10203 q = u*h*w = 20m^3/s10201 w = 3 m (width of channel) 10202 10203 q = u*h*w = 12 m^3/s 10204 10204 10205 10205 … … 10243 10243 10244 10244 Br = Reflective_boundary(domain) # Side walls 10245 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 5m inlet:10245 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 10246 10246 10247 10247 … … 10300 10300 u = 2 m/s 10301 10301 h = 1 m 10302 w = 5m (width of channel)10303 10304 q = u*h*w = 10m^3/s10302 w = 3 m (width of channel) 10303 10304 q = u*h*w = 6 m^3/s 10305 10305 Es = h + 0.5*v*v/g # Specific energy head [m] 10306 10306 Et = w + 0.5*v*v/g # Total energy head [m] … … 10346 10346 10347 10347 Br = Reflective_boundary(domain) # Side walls 10348 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 5m inlet:10348 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 10349 10349 10350 10350 -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r5450 r5729 1406 1406 1407 1407 1408 1409 def test_get_flow_through_cross_section_with_geo(self): 1410 """test_get_flow_through_cross_section(self): 1411 1412 Test that the total flow through a cross section can be 1413 correctly obtained at run-time from the ANUGA domain. 1414 1415 This test creates a flat bed with a known flow through it and tests 1416 that the function correctly returns the expected flow. 1417 1418 The specifics are 1419 e = -1 m 1420 u = 2 m/s 1421 h = 2 m 1422 w = 3 m (width of channel) 1423 1424 q = u*h*w = 12 m^3/s 1425 1426 1427 This run tries it with georeferencing and with elevation = -1 1428 1429 """ 1430 1431 import time, os 1432 from Numeric import array, zeros, allclose, Float, concatenate 1433 from Scientific.IO.NetCDF import NetCDFFile 1434 1435 # Setup 1436 from mesh_factory import rectangular 1437 1438 # Create basic mesh (20m x 3m) 1439 width = 3 1440 length = 20 1441 t_end = 1 1442 points, vertices, boundary = rectangular(length, width, 1443 length, width) 1444 1445 # Create shallow water domain 1446 domain = Domain(points, vertices, boundary, 1447 geo_reference=Geo_reference(56,308500,6189000)) 1448 1449 domain.default_order = 2 1450 domain.set_quantities_to_be_stored(None) 1451 1452 1453 e = -1.0 1454 w = 1.0 1455 h = w-e 1456 u = 2.0 1457 uh = u*h 1458 1459 Br = Reflective_boundary(domain) # Side walls 1460 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1461 1462 1463 # Initial conditions 1464 domain.set_quantity('elevation', e) 1465 domain.set_quantity('stage', w) 1466 domain.set_quantity('xmomentum', uh) 1467 domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1468 1469 1470 # Interpolation points down the middle 1471 I = [[0, width/2.], 1472 [length/2., width/2.], 1473 [length, width/2.]] 1474 interpolation_points = domain.geo_reference.get_absolute(I) 1475 1476 # Shortcuts to quantites 1477 stage = domain.get_quantity('stage') 1478 xmomentum = domain.get_quantity('xmomentum') 1479 ymomentum = domain.get_quantity('ymomentum') 1480 1481 for t in domain.evolve(yieldstep=0.1, finaltime = t_end): 1482 # Check that quantities are they should be in the interior 1483 1484 w_t = stage.get_values(interpolation_points) 1485 uh_t = xmomentum.get_values(interpolation_points) 1486 vh_t = ymomentum.get_values(interpolation_points) 1487 1488 assert allclose(w_t, w) 1489 assert allclose(uh_t, uh) 1490 assert allclose(vh_t, 0.0) 1491 1492 1493 # Check flows through the middle 1494 for i in range(5): 1495 x = length/2. + i*0.23674563 # Arbitrary 1496 cross_section = [[x, 0], [x, width]] 1497 1498 cross_section = domain.geo_reference.get_absolute(cross_section) 1499 Q = domain.get_flow_through_cross_section(cross_section, 1500 verbose=False) 1501 1502 assert allclose(Q, uh*width) 1503 1504 1505 1506 1507 1408 1508 1409 1509 def test_another_runup_example(self): … … 1431 1531 domain = Domain(points, vertices, boundary) # Create domain 1432 1532 domain.set_quantities_to_be_stored(None) 1433 domain.set_maximum_allowed_speed(100) # 1533 domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this 1434 1534 1435 1535 # FIXME (Ole): Need tests where this is commented out 1436 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 1536 domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7) 1437 1537 domain.H0 = 0 # Backwards compatibility (6/2/7) 1438 1538 domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8) … … 5805 5905 if __name__ == "__main__": 5806 5906 5807 suite = unittest.makeSuite(Test_Shallow_Water,'test')5808 5809 #suite = unittest.makeSuite(Test_Shallow_Water,'test_bedslope_problem_first_order_moresteps')5907 #suite = unittest.makeSuite(Test_Shallow_Water,'test') 5908 5909 suite = unittest.makeSuite(Test_Shallow_Water,'test_get_flow_through_cross_section_with_geo') 5810 5910 #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain') 5811 5911 #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters') -
anuga_core/source/anuga/utilities/numerical_tools.py
r5561 r5729 313 313 314 314 return bins 315 316 315 317 316 318 def get_machine_precision():
Note: See TracChangeset
for help on using the changeset viewer.