Changeset 664

Timestamp:

Dec 2, 2004, 3:28:33 PM (20 years ago)

Author:

ole

Message:

Added arbitrary functions to Polygon_function

Location:

inundation/ga/storm_surge/pyvolution

Files:

• flatbed.py (modified) (4 diffs)
• test_util.py (modified) (5 diffs)
• util.py (modified) (7 diffs)

inundation/ga/storm_surge/pyvolution/flatbed.py

@@ -17 +17 @@
 
 #Create basic mesh
-N = 20
+N = 50
 points, vertices, boundary = rectangular(N, N, 100, 100)
 
@@ -25 +25 @@
 domain.set_name('polygons')
 print "Output being written to " + data_dir + sep + \
-              domain.filename + "." + domain.format
+              domain.filename + "_size%d." %len(domain) + domain.format
         
 
@@ -36 +36 @@
 domain.set_quantity('friction', manning)
 
+def wiggle(x, y):
+    from Numeric import sin
+    from math import pi
+    return 10 + sin(2*pi*x/10)
 
-
+def slope(x, y):
+    return 20*(x/100+y/100)    
 
 #Define polynomials
@@ -47 +52 @@
 #Set elevation  
 domain.set_quantity('elevation', 
-        Polygon_function([(p0, 23), (p1,10), (p2,15)]))
+        Polygon_function([(p0,slope), (p1,wiggle), (p2,15)]))
+#domain.set_quantity('level', 
+#       Polygon_function([(p0,slope), (p1,wiggle), (p2,15)]))   
 
 

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -8 +8 @@
 from util import *
 from config import epsilon
+
+
+def test_function(x, y):
+    return x+y
 
 class TestCase(unittest.TestCase):
@@ -334 +338 @@
         
         
-    def test_polygon_function(self):
+    def test_polygon_function_constants(self):
         p1 = [[0,0], [10,0], [10,10], [0,10]]
         p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
@@ -353 +357 @@
         assert allclose(z, [2,1,0,2])   
         
+        
+    def test_polygon_function_callable(self):
+        """Check that values passed into Polygon_function can be callable 
+        themselves.
+        """
+        p1 = [[0,0], [10,0], [10,10], [0,10]]
+        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
+        
+        f = Polygon_function( [(p1, test_function)] )   
+        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
+        assert allclose(z, [10,14,0,0])
+                
+        #Combined       
+        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
+        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
+        assert allclose(z, [2,14,0,2])  
+        
+        
+        #Combined w default     
+        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
+        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
+        assert allclose(z, [2,14,3.14,2])               
+        
+        
+        #Combined w default func        
+        f = Polygon_function( [(p1, test_function), (p2, 2.0)], 
+                              default = test_function)
+        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
+        assert allclose(z, [2,14,35,2])                 
         
         
@@ -467 +500 @@
         assert inside_polygon( (10.5, 10.5), polygon )                  
         
-        #FIXME: Fails
-        #assert not inside_polygon( (5.5, 5.5), polygon )                               
+        #FIXME: Fails if point is 5.5, 5.5
+        assert not inside_polygon( (0, 5.5), polygon )                          
 
         #Polygon with a hole    
@@ -475 +508 @@
                         
         assert inside_polygon( (0, -0.5), polygon )
-        
-        #FIXME: Fails   
-        #assert not inside_polygon( (0.5, 0.5), polygon )       
+        assert not inside_polygon( (0.5, 0.5), polygon )        
         
 

inundation/ga/storm_surge/pyvolution/util.py

@@ -57 +57 @@
     """
     
-    #FIXME: Could do with some optimisation         
+    #FIXME: Could do with some C-optimisation         
         
     from Numeric import array, dot, allclose
@@ -94 +94 @@
        point - Tuple of (x, y) coordinates, or list of tuples
        polygon - list of vertices of polygon
+       closed - (optional) determine whether points on boundary should be 
+       regarded as belonging to the polygon (closed = True) 
+       or not (closed = False) 
     
     Output:
@@ -109 +112 @@
        
     Remarks:
-       <Copy from source Franklins code >    
+       The vertices may be listed clockwise or counterclockwise and
+       the first point may optionally be repeated.    
+       Polygons do not need to be convex.
+       Polygons can have holes in them and points inside a hole is 
+       regarded as being outside the polygon.          
                
-               
-       points at a vertex are considered to be inside          
        
-    Algorithm due to Darel Finley, http://www.alienryderflex.com/polygon/ 
-    
-    
-    FIXME: More comments about algoritm
-    FIXME: Deal with case where point on border of polygon (closed)
-       
+    Algorithm is based on work by Darel Finley, 
+    http://www.alienryderflex.com/polygon/ 
+    
     """    
 
@@ -166 +168 @@
         y = points[k, 1]        
 
-        
         inside = False
         for i in range(N):
@@ -187 +188 @@
         if inside: indices.append(k)
 
-        #print x, y, polygon, inside
-
-    
     if one_point: 
         return inside
@@ -222 +220 @@
     Note: If two polygons overlap, the one last in the list takes precedence
 
-    Note: default can only be a constant.    
     """
     
@@ -253 +250 @@
         x = array(x).astype(Float)
         y = array(y).astype(Float)      
-        x = reshape(x, (len(x), 1))
-        y = reshape(y, (len(y), 1))     
-                
-        points = concatenate( (x,y), axis=1 )
-        
-        z = ones(x.shape[0], Float) * self.default
+        
+        N = len(x)
+        assert len(y) == N
+        
+        points = concatenate( (reshape(x, (N, 1)), 
+                               reshape(y, (N, 1))), axis=1 )
+
+        if callable(self.default):
+            z = self.default(x,y)
+        else:                   
+            z = ones(N, Float) * self.default
+            
         for polygon, value in self.regions:
             indices = inside_polygon(points, polygon)
-            for i in indices:
-                z[i] = value 
-
-        #from sys import exit; exit()   
+            
+            #FIXME: This needs to be vectorised 
+            if callable(value):
+                for i in indices:
+                    xx = array([x[i]])
+                    yy = array([y[i]])
+                    z[i] = value(xx, yy)[0]
+            else:    
+                for i in indices:
+                    z[i] = value 
+
         return z