Changeset 5570 for anuga_core/source

Timestamp:

Jul 25, 2008, 1:35:45 PM (17 years ago)

Author:

ole

Message:

Added input tests regarding polygons and points

Location:

anuga_core/source/anuga

Files:

• culvert_flows/culvert_class.py (modified) (6 diffs)
• shallow_water/shallow_water_domain.py (modified) (4 diffs)
• utilities/polygon.py (modified) (1 diff)
• utilities/polygon_ext.c (modified) (1 diff)

anuga_core/source/anuga/culvert_flows/culvert_class.py

@@ -2 +2 @@
 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
 from anuga.utilities.system_tools import log_to_file
+from anuga.utilities.polygon import inside_polygon
+from anuga.utilities.polygon import is_inside_polygon
+
 
 class Culvert_flow:
@@ -125 +128 @@
             #               P['enquiry_polygon1']],
             #              figname='culvert_polygon_output')
-        
+
+
+        # Check that all polygons lie within the mesh.
+        bounding_polygon = domain.get_boundary_polygon()
+        for key in P.keys():
+            print 'Key', key
+            for point in P[key]:
+
+                print 'Passing in:', point
+                msg = 'Point %s in polygon %s for culvert %s did not'\
+                      %(str(point), key, self.label)
+                msg += 'fall within the domain boundary.'
+                assert is_inside_polygon(point, bounding_polygon), msg
+        
+
+        # Create inflow object at each end of the culvert. 
         self.openings = []
         self.openings.append(Inflow(domain,
@@ -188 +206 @@
     def __call__(self, domain):
         from anuga.utilities.numerical_tools import mean
-        from anuga.utilities.polygon import inside_polygon
+        
         from anuga.config import g, epsilon
         from Numeric import take, sqrt
@@ -209 +227 @@
         openings = self.openings   # There are two Opening [0] and [1]
         for i, opening in enumerate(openings):
-            stage = domain.quantities['stage'].get_values(location='centroids',
-                                                          indices=opening.exchange_indices)
-            elevation = domain.quantities['elevation'].get_values(location='centroids',
-                                                                  indices=opening.exchange_indices)
+            dq = domain.quantities
+            
+            stage = dq['stage'].get_values(location='centroids',
+                                           indices=opening.exchange_indices)
+            elevation = dq['elevation'].get_values(location='centroids',
+                                                   indices=opening.exchange_indices)
 
             # Indices corresponding to energy enquiry field for this opening
             coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)
-            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 
+            enquiry_indices = inside_polygon(coordinates,
+                                             self.enquiry_polygons[i]) 
 
             if len(enquiry_indices) == 0:
@@ -223 +244 @@
             
             # Get model values for points in enquiry polygon for this opening
-            dq = domain.quantities
-            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
-            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
-            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                        
-            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
+            stage = dq['stage'].get_values(location='centroids',
+                                           indices=enquiry_indices)
+            xmomentum = dq['xmomentum'].get_values(location='centroids',
+                                                   indices=enquiry_indices)
+            ymomentum = dq['ymomentum'].get_values(location='centroids',
+                                                   indices=enquiry_indices)
+            elevation = dq['elevation'].get_values(location='centroids',
+                                                   indices=enquiry_indices)
             depth = stage - elevation
 
@@ -234 +258 @@
             ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
             uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
+            print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum
             v = mean(sqrt(ux**2+uy**2))      # Average velocity
             w = mean(stage)                  # Average stage

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -108 +108 @@
 
 
-from anuga.utilities.polygon import inside_polygon, polygon_area        
+from anuga.utilities.polygon import inside_polygon, polygon_area, is_inside_polygon
 
 
@@ -1501 +1501 @@
             
                      
-        from math import pi
+        from math import pi, cos, sin
 
         self.domain = domain
@@ -1513 +1513 @@
                          # previous timestep in order to obtain rate
 
+                         
+        bounding_polygon = domain.get_boundary_polygon()
+
+
         # Update area if applicable
         self.exchange_area = None        
@@ -1521 +1525 @@
 
             self.exchange_area = radius**2*pi
+
+            # Check that circle center lies within the mesh.
+            msg = 'Center %s specified for forcing term did not' %(str(center))
+            msg += 'fall within the domain boundary.'
+            assert is_inside_polygon(center, bounding_polygon), msg
+
+            # Check that circle periphery lies within the mesh.
+            N = 100
+            periphery_points = []
+            for i in range(N):
+
+                theta = 2*pi*i/100
+                
+                x = center[0] + radius*cos(theta)
+                y = center[1] + radius*sin(theta)
+
+                periphery_points.append([x,y])
+
+
+            for point in periphery_points:
+                msg = 'Point %s on periphery for forcing term did not' %(str(point))
+                msg += ' fall within the domain boundary.'
+                assert is_inside_polygon(point, bounding_polygon), msg
+
         
         if polygon is not None:
             self.exchange_area = polygon_area(self.polygon)
+
+            # Check that polygon lies within the mesh.
+            for point in self.polygon:
+                msg = 'Point %s in polygon for forcing term did not' %(str(point))
+                msg += 'fall within the domain boundary.'
+                assert is_inside_polygon(point, bounding_polygon), msg
+        
+
+
+            
 
 

anuga_core/source/anuga/utilities/polygon.py

@@ -457 +457 @@
         raise msg
 
-    #if verbose: print 'check'
-
-    assert len(polygon.shape) == 2,\
-       'Polygon array must be a 2d array of vertices'
-
-    assert polygon.shape[1] == 2,\
-       'Polygon array must have two columns'
-
-    assert len(points.shape) == 2,\
-       'Points array must be a 2d array'
-
-    assert points.shape[1] == 2,\
-       'Points array must have two columns'
+    msg = 'Polygon array must be a 2d array of vertices'
+    assert len(polygon.shape) == 2, msg
+
+    msg = 'Polygon array must have two columns' 
+    assert polygon.shape[1] == 2, msg
+
+
+    msg = 'Points array must be 1 or 2 dimensional.'
+    msg += ' I got %d dimensions' %len(points.shape)
+    assert 0 < len(points.shape) < 3, msg
+
+
+    if len(points.shape) == 1:
+        # Only one point was passed in.
+        # Convert to array of points
+        points = reshape(points, (1,2))
+
+    
+    msg = 'Point array must have two columns (x,y), '
+    msg += 'I got points.shape[1] == %d' %points.shape[0]
+    assert points.shape[1] == 2, msg
+
+       
+    msg = 'Points array must be a 2d array. I got %s' %str(points[:30])
+    assert len(points.shape) == 2, msg
+
+    msg = 'Points array must have two columns'
+    assert points.shape[1] == 2, msg
+
 
     N = polygon.shape[0] #Number of vertices in polygon

anuga_core/source/anuga/utilities/polygon_ext.c

@@ -199 +199 @@
   //Find min and max of poly used for optimisation when points
   //are far away from polygon
+  
+  //FIXME(Ole): Pass in rtol and atol from Python
 
   minpx = polygon[0]; maxpx = minpx;