Changeset 5570
 Timestamp:
 Jul 25, 2008, 1:35:45 PM (11 years ago)
 Location:
 anuga_core/source/anuga
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/culvert_flows/culvert_class.py
r5453 r5570 2 2 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons 3 3 from anuga.utilities.system_tools import log_to_file 4 from anuga.utilities.polygon import inside_polygon 5 from anuga.utilities.polygon import is_inside_polygon 6 4 7 5 8 class Culvert_flow: … … 125 128 # P['enquiry_polygon1']], 126 129 # figname='culvert_polygon_output') 127 130 131 132 # Check that all polygons lie within the mesh. 133 bounding_polygon = domain.get_boundary_polygon() 134 for key in P.keys(): 135 print 'Key', key 136 for point in P[key]: 137 138 print 'Passing in:', point 139 msg = 'Point %s in polygon %s for culvert %s did not'\ 140 %(str(point), key, self.label) 141 msg += 'fall within the domain boundary.' 142 assert is_inside_polygon(point, bounding_polygon), msg 143 144 145 # Create inflow object at each end of the culvert. 128 146 self.openings = [] 129 147 self.openings.append(Inflow(domain, … … 188 206 def __call__(self, domain): 189 207 from anuga.utilities.numerical_tools import mean 190 from anuga.utilities.polygon import inside_polygon208 191 209 from anuga.config import g, epsilon 192 210 from Numeric import take, sqrt … … 209 227 openings = self.openings # There are two Opening [0] and [1] 210 228 for i, opening in enumerate(openings): 211 stage = domain.quantities['stage'].get_values(location='centroids', 212 indices=opening.exchange_indices) 213 elevation = domain.quantities['elevation'].get_values(location='centroids', 214 indices=opening.exchange_indices) 229 dq = domain.quantities 230 231 stage = dq['stage'].get_values(location='centroids', 232 indices=opening.exchange_indices) 233 elevation = dq['elevation'].get_values(location='centroids', 234 indices=opening.exchange_indices) 215 235 216 236 # Indices corresponding to energy enquiry field for this opening 217 237 coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y) 218 enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 238 enquiry_indices = inside_polygon(coordinates, 239 self.enquiry_polygons[i]) 219 240 220 241 if len(enquiry_indices) == 0: … … 223 244 224 245 # Get model values for points in enquiry polygon for this opening 225 dq = domain.quantities 226 stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices) 227 xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices) 228 ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices) 229 elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices) 246 stage = dq['stage'].get_values(location='centroids', 247 indices=enquiry_indices) 248 xmomentum = dq['xmomentum'].get_values(location='centroids', 249 indices=enquiry_indices) 250 ymomentum = dq['ymomentum'].get_values(location='centroids', 251 indices=enquiry_indices) 252 elevation = dq['elevation'].get_values(location='centroids', 253 indices=enquiry_indices) 230 254 depth = stage  elevation 231 255 … … 234 258 ux = xmomentum/(depth+velocity_protection/depth) # Velocity (xdirection) 235 259 uy = ymomentum/(depth+velocity_protection/depth) # Velocity (ydirection) 260 print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum 236 261 v = mean(sqrt(ux**2+uy**2)) # Average velocity 237 262 w = mean(stage) # Average stage 
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r5506 r5570 108 108 109 109 110 from anuga.utilities.polygon import inside_polygon, polygon_area 110 from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon 111 111 112 112 … … 1501 1501 1502 1502 1503 from math import pi 1503 from math import pi, cos, sin 1504 1504 1505 1505 self.domain = domain … … 1513 1513 # previous timestep in order to obtain rate 1514 1514 1515 1516 bounding_polygon = domain.get_boundary_polygon() 1517 1518 1515 1519 # Update area if applicable 1516 1520 self.exchange_area = None … … 1521 1525 1522 1526 self.exchange_area = radius**2*pi 1527 1528 # Check that circle center lies within the mesh. 1529 msg = 'Center %s specified for forcing term did not' %(str(center)) 1530 msg += 'fall within the domain boundary.' 1531 assert is_inside_polygon(center, bounding_polygon), msg 1532 1533 # Check that circle periphery lies within the mesh. 1534 N = 100 1535 periphery_points = [] 1536 for i in range(N): 1537 1538 theta = 2*pi*i/100 1539 1540 x = center[0] + radius*cos(theta) 1541 y = center[1] + radius*sin(theta) 1542 1543 periphery_points.append([x,y]) 1544 1545 1546 for point in periphery_points: 1547 msg = 'Point %s on periphery for forcing term did not' %(str(point)) 1548 msg += ' fall within the domain boundary.' 1549 assert is_inside_polygon(point, bounding_polygon), msg 1550 1523 1551 1524 1552 if polygon is not None: 1525 1553 self.exchange_area = polygon_area(self.polygon) 1554 1555 # Check that polygon lies within the mesh. 1556 for point in self.polygon: 1557 msg = 'Point %s in polygon for forcing term did not' %(str(point)) 1558 msg += 'fall within the domain boundary.' 1559 assert is_inside_polygon(point, bounding_polygon), msg 1560 1561 1562 1563 1526 1564 1527 1565 
anuga_core/source/anuga/utilities/polygon.py
r5403 r5570 457 457 raise msg 458 458 459 #if verbose: print 'check' 460 461 assert len(polygon.shape) == 2,\ 462 'Polygon array must be a 2d array of vertices' 463 464 assert polygon.shape[1] == 2,\ 465 'Polygon array must have two columns' 466 467 assert len(points.shape) == 2,\ 468 'Points array must be a 2d array' 469 470 assert points.shape[1] == 2,\ 471 'Points array must have two columns' 459 msg = 'Polygon array must be a 2d array of vertices' 460 assert len(polygon.shape) == 2, msg 461 462 msg = 'Polygon array must have two columns' 463 assert polygon.shape[1] == 2, msg 464 465 466 msg = 'Points array must be 1 or 2 dimensional.' 467 msg += ' I got %d dimensions' %len(points.shape) 468 assert 0 < len(points.shape) < 3, msg 469 470 471 if len(points.shape) == 1: 472 # Only one point was passed in. 473 # Convert to array of points 474 points = reshape(points, (1,2)) 475 476 477 msg = 'Point array must have two columns (x,y), ' 478 msg += 'I got points.shape[1] == %d' %points.shape[0] 479 assert points.shape[1] == 2, msg 480 481 482 msg = 'Points array must be a 2d array. I got %s' %str(points[:30]) 483 assert len(points.shape) == 2, msg 484 485 msg = 'Points array must have two columns' 486 assert points.shape[1] == 2, msg 487 472 488 473 489 N = polygon.shape[0] #Number of vertices in polygon 
anuga_core/source/anuga/utilities/polygon_ext.c
r5403 r5570 199 199 //Find min and max of poly used for optimisation when points 200 200 //are far away from polygon 201 202 //FIXME(Ole): Pass in rtol and atol from Python 201 203 202 204 minpx = polygon[0]; maxpx = minpx;
Note: See TracChangeset
for help on using the changeset viewer.