Index: /inundation/ga/storm_surge/pyvolution/flatbed.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/flatbed.py	(revision 663)
+++ /inundation/ga/storm_surge/pyvolution/flatbed.py	(revision 664)
@@ -17,5 +17,5 @@
 
 #Create basic mesh
-N = 20
+N = 50
 points, vertices, boundary = rectangular(N, N, 100, 100)
 
@@ -25,5 +25,5 @@
 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,6 +36,11 @@
 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,5 +52,7 @@
 #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)]))	
 
 
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 663)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 664)
@@ -8,4 +8,8 @@
 from util import *
 from config import epsilon
+
+
+def test_function(x, y):
+    return x+y
 
 class TestCase(unittest.TestCase):
@@ -334,5 +338,5 @@
 	
 	
-    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,4 +357,33 @@
 	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,6 +500,6 @@
         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,7 +508,5 @@
 		   	
         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 )	
 	
 
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 663)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 664)
@@ -57,5 +57,5 @@
     """
     
-    #FIXME: Could do with some optimisation         
+    #FIXME: Could do with some C-optimisation         
 	
     from Numeric import array, dot, allclose
@@ -94,4 +94,7 @@
        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,15 +112,14 @@
        
     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,5 +168,4 @@
 	y = points[k, 1]	
 
-	
         inside = False
         for i in range(N):
@@ -187,7 +188,4 @@
         if inside: indices.append(k)
 
-	#print x, y, polygon, inside
-
-    
     if one_point: 
         return inside
@@ -222,5 +220,4 @@
     Note: If two polygons overlap, the one last in the list takes precedence
 
-    Note: default can only be a constant.    
     """
     
@@ -253,16 +250,29 @@
 	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    		 
 
