Changeset 5570


Ignore:
Timestamp:
Jul 25, 2008, 1:35:45 PM (11 years ago)
Author:
ole
Message:

Added input tests regarding polygons and points

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r5453 r5570  
    22from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
    33from anuga.utilities.system_tools import log_to_file
     4from anuga.utilities.polygon import inside_polygon
     5from anuga.utilities.polygon import is_inside_polygon
     6
    47
    58class Culvert_flow:
     
    125128            #               P['enquiry_polygon1']],
    126129            #              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.
    128146        self.openings = []
    129147        self.openings.append(Inflow(domain,
     
    188206    def __call__(self, domain):
    189207        from anuga.utilities.numerical_tools import mean
    190         from anuga.utilities.polygon import inside_polygon
     208       
    191209        from anuga.config import g, epsilon
    192210        from Numeric import take, sqrt
     
    209227        openings = self.openings   # There are two Opening [0] and [1]
    210228        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)
    215235
    216236            # Indices corresponding to energy enquiry field for this opening
    217237            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])
    219240
    220241            if len(enquiry_indices) == 0:
     
    223244           
    224245            # 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)
    230254            depth = stage - elevation
    231255
     
    234258            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
    235259            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
     260            print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
    236261            v = mean(sqrt(ux**2+uy**2))      # Average velocity
    237262            w = mean(stage)                  # Average stage
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5506 r5570  
    108108
    109109
    110 from anuga.utilities.polygon import inside_polygon, polygon_area       
     110from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
    111111
    112112
     
    15011501           
    15021502                     
    1503         from math import pi
     1503        from math import pi, cos, sin
    15041504
    15051505        self.domain = domain
     
    15131513                         # previous timestep in order to obtain rate
    15141514
     1515                         
     1516        bounding_polygon = domain.get_boundary_polygon()
     1517
     1518
    15151519        # Update area if applicable
    15161520        self.exchange_area = None       
     
    15211525
    15221526            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
    15231551       
    15241552        if polygon is not None:
    15251553            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           
    15261564
    15271565
  • anuga_core/source/anuga/utilities/polygon.py

    r5403 r5570  
    457457        raise msg
    458458
    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
    472488
    473489    N = polygon.shape[0] #Number of vertices in polygon
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r5403 r5570  
    199199  //Find min and max of poly used for optimisation when points
    200200  //are far away from polygon
     201 
     202  //FIXME(Ole): Pass in rtol and atol from Python
    201203
    202204  minpx = polygon[0]; maxpx = minpx;
Note: See TracChangeset for help on using the changeset viewer.