Index: /inundation/ga/storm_surge/examples/beach.py
===================================================================
--- /inundation/ga/storm_surge/examples/beach.py	(revision 658)
+++ /inundation/ga/storm_surge/examples/beach.py	(revision 659)
@@ -21,4 +21,5 @@
 from pmesh2domain import pmesh_to_domain_instance
 from config import data_dir
+from util import read_polygon, Polygon_function
 
 ######################
@@ -49,41 +50,13 @@
 
 
-def inside(p,polygon):
-    return True
 
-def friction(x, y):
-    """Read friction from polygon file and set accordingly
-    """
 
-    from Numeric import array
-    x = array(x).astype('f')
-    y = array(y).astype('f')    
+p0 = [[20,20], [40,20], [40,40], [20,40]]         
     
-    #Get polygon
-    fid = open('beach.poly')
-    lines = fid.readlines()
-    fid.close()
-    polygon = []
-    for line in lines:
-        fields = line.split(',')
-        polygon.append( [float(fields[0]), float(fields[1])] )
-
-    print polygon
-
-    z = 0*x
-    for i in range(len(x)):
-        if inside( (x[i], y[i]), polygon):
-            z[i] = 1
-        else:
-            z[i] = 0
-
-friction([0],[1])    
-
-domain.set_quantity('elevation', x_slope)
-domain.set_quantity('level', x_slope)
+domain.set_quantity('elevation', Polygon_function( [(p0, 50)]), 'centroids' )
+domain.set_quantity('level', Polygon_function( [(p0, 50)] ), 'centroids' )
 domain.set_quantity('friction', 0.07)
 
 print domain.get_extent()
-import sys; sys.exit()
 
 ######################
Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 659)
@@ -121,4 +121,5 @@
             self.set_quantity(key, quantity_dict[key], location='vertices')
         
+
     def set_quantity(self, name, X, location='vertices', indexes = None):
         """Set values for named quantity
Index: /inundation/ga/storm_surge/pyvolution/flatbed.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/flatbed.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/flatbed.py	(revision 659)
@@ -10,17 +10,17 @@
 from mesh_factory import rectangular
 from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
-     Constant_height, wind_Stress
+     Constant_height
 from Numeric import array
+from util import Polygon_function, read_polygon
 
 #Create basic mesh
-N = 150
-points, vertices, boundary = rectangular(N, N, 1000, 1000)
+N = 20
+points, vertices, boundary = rectangular(N, N, 100, 100)
 
 #Create shallow water domain
 domain = Domain(points, vertices, boundary)
-domain.smooth = True
-domain.visualise = False
 domain.store = True
-domain.default_order=1
+domain.set_name('polygons')
+domain.default_order=2
 
 #Set driving forces
@@ -29,4 +29,20 @@
 inflow_level = 10.0
 domain.set_quantity('friction', manning)
+
+
+
+
+#Define polynomials
+p0 = [[20,27], [30,25], [40,40], [20,40]]         
+p1 = [[80,19], [90,20], [85,50], [80,55], [75,58], [70,60], [60,24]]         
+p2 = read_polygon('testpoly.txt')
+
+
+#Set elevation  
+domain.set_quantity('elevation', 
+	Polygon_function([(p0, 23), (p1,10), (p2,15)]))
+
+
+ 	  
 
 
@@ -43,5 +59,5 @@
 ######################
 #Evolution
-for t in domain.evolve(yieldstep = 10, finaltime = 500):
+for t in domain.evolve(yieldstep = 1, finaltime = 1):
     domain.write_time()
 
Index: /inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 659)
@@ -147,5 +147,8 @@
             #Intialise centroid and edge_values
             self.interpolate()
-            
+	    
+        if location == 'centroids':
+            #Extrapolate 1st order - to capture notion of area being specified
+            self.extrapolate_first_order()            
           
     def get_values(self, location='vertices', indexes = None):
@@ -479,4 +482,15 @@
 
             
+    def extrapolate_first_order(self):
+        """Extrapolate conserved quantities from centroid to
+        vertices for each volume using
+        first order scheme.
+        """
+        
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        for i in range(3):
+            qv[:,i] = qc
             
 
@@ -486,5 +500,5 @@
 
     boundary values, storage and method for updating, and
-    methods for extrapolation from centropid to vertices inluding
+    methods for (second order) extrapolation from centroid to vertices inluding
     gradients and limiters
     """
@@ -525,17 +539,4 @@
         
         
-    def extrapolate_first_order(self):
-        """Extrapolate conserved quantities from centroid to
-        vertices for each volume using
-        first order scheme.
-        """
-        
-        qc = self.centroid_values
-        qv = self.vertex_values
-
-        for i in range(3):
-            qv[:,i] = qc
-            
-            
     def extrapolate_second_order(self):
         #Call correct module function
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 659)
@@ -81,13 +81,13 @@
         bed = self.quantities['elevation']        
 
-        msg = 'All water levels must be greater than the bed elevation'
-        assert alltrue( greater_equal(
-            level.vertex_values, bed.vertex_values )), msg
-
-        assert alltrue( greater_equal(
-            level.edge_values, bed.edge_values )), msg
-
-        assert alltrue( greater_equal(
-            level.centroid_values, bed.centroid_values )), msg        
+        #msg = 'All water levels must be greater than the bed elevation'
+        #assert alltrue( greater_equal(
+        #    level.vertex_values, bed.vertex_values )), msg
+        # 
+        #assert alltrue( greater_equal(
+        #    level.edge_values, bed.edge_values )), msg
+        # 
+        #assert alltrue( greater_equal(
+        #    level.centroid_values, bed.centroid_values )), msg        
 
 
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 659)
@@ -356,5 +356,5 @@
 	
     def test_point_on_line(self):
-        	
+            	
 	#Endpoints first	
 	assert point_on_line( (0, 0), [(0,0), (1,0)] ) 	
@@ -371,8 +371,17 @@
 	assert not point_on_line( (0, 0.5), [(0,0), (1,1)] ) 
 	
+	#From real example that failed
+	assert not point_on_line( (40,50), [(40,20), (40,40)] )
+	
+	
+	#From real example that failed
+	assert not point_on_line( (40,19), [(40,20), (40,40)] )	
+	 		
+	
 	
 	
     def test_inside_polygon(self):
-    
+
+        
         #Simplest case: Polygon is the unit square
         polygon = [[0,0], [1,0], [1,1], [0,1]]
@@ -394,5 +403,20 @@
 	assert not inside_polygon( (0.5, 0.), polygon, closed=False)	
 	assert not inside_polygon( (1., 0.5), polygon, closed=False)	
-	
+
+
+    
+        #From real example (that failed)
+	polygon = [[20,20], [40,20], [40,40], [20,40]]         
+	points = [ [40, 50] ]
+	res = inside_polygon(points, polygon)
+	assert len(res) == 0
+	
+	polygon = [[20,20], [40,20], [40,40], [20,40]]         
+        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
+	res = inside_polygon(points, polygon)
+	assert len(res) == 2
+	assert allclose(res, [0,1])
+        
+		
 	      
 	#More convoluted and non convex polygon
@@ -455,6 +479,4 @@
         #assert not inside_polygon( (0.5, 0.5), polygon )	
 	
-	points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]			
-	
 
 		      
Index: /inundation/ga/storm_surge/pyvolution/testpoly.txt
===================================================================
--- /inundation/ga/storm_surge/pyvolution/testpoly.txt	(revision 659)
+++ /inundation/ga/storm_surge/pyvolution/testpoly.txt	(revision 659)
@@ -0,0 +1,7 @@
+40, 70
+50, 70
+60, 80
+50, 90
+45, 95
+40, 90
+30, 80
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 658)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 659)
@@ -70,19 +70,22 @@
     y1 = line[1][1]
 
-    
-    a = array([x - x0, y - y0])        		
+    a = array([x - x0, y - y0]) 
+    a_normal = array([a[1], -a[0]])
+           		
     b = array([x1 - x0, y1 - y0])
     
-    len_a = sqrt(sum(a**2))    	       
-    len_b = sqrt(sum(b**2))    	       	    
-        	
-    if len_a == 0 or len_a == len_b:
-       return True
+    if dot(a_normal, b) == 0:
+        #Point is somewhere on the infinite extension of the line
+
+        len_a = sqrt(sum(a**2))    	       
+        len_b = sqrt(sum(b**2))    	       	                
+        if dot(a, b) >= 0 and len_a <= len_b:
+           return True
+	else:   
+           return False
     else:
-       x = dot(a, b)/len_b       
-       return allclose(x, len_a)
-
-        
-    
+      return False  
+    
+
 def inside_polygon(point, polygon, closed = True):
     """Determine whether points are inside or outside a polygon
@@ -155,5 +158,5 @@
     py = polygon[:,1]	
 
-        	        
+    
     #Begin algorithm
     indices = []
@@ -162,4 +165,5 @@
 	x = points[k, 0]
 	y = points[k, 1]	
+
 	
         inside = False
@@ -257,9 +261,26 @@
 	for polygon, value in self.regions:
 	    indices = inside_polygon(points, polygon)
-	    
 	    for i in indices:
                 z[i] = value 
 
+	#from sys import exit; exit()	
         return z    		 
+
+def read_polygon(filename):
+    """Read points assumed to form a polygon
+       There must be exactly two numbers in each line
+    """		     
+    
+    #Get polygon
+    fid = open(filename)
+    lines = fid.readlines()
+    fid.close()
+    polygon = []
+    for line in lines:
+        fields = line.split(',')
+        polygon.append( [float(fields[0]), float(fields[1])] )
+
+    return polygon	
+	     
 	     
 class File_function:
