Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py	(revision 4554)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py	(revision 4555)
@@ -141,5 +141,6 @@
     #                            rtol=rtol, atol=atol)
     if plotting is False:
-        assert allclose(reference_value, computed_value, rtol=rtol, atol=atol), msg
+        assert allclose(reference_value, computed_value,
+                        rtol=rtol, atol=atol), msg
     
 
@@ -215,2 +216,28 @@
 
 
+# Check max runup
+
+from anuga.shallow_water.data_manager import get_maximum_inundation_elevation
+from anuga.shallow_water.data_manager import get_maximum_inundation_location
+from anuga.utilities.polygon import is_inside_polygon
+
+q = get_maximum_inundation_elevation(sww_filename)
+loc = get_maximum_inundation_location(sww_filename)
+
+
+print 'Max runup elevation: ', q
+print 'Max runup location:  ', loc
+
+
+from create_okushiri import gulleys
+assert is_inside_polygon(loc, gulleys)
+
+# FIXME more asserts here
+
+
+
+#msg = 'We got %f, should have been %f' %(q, q_max)
+#assert allclose(q, q_max, rtol=1.0/N), msg
+##print 'loc', loc, q
+#assert allclose(-loc[0]/2, q) # From topography formula 
+
Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py	(revision 4554)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py	(revision 4555)
@@ -11,4 +11,80 @@
 
 import project
+
+
+
+
+#--------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with a
+# tagged
+# boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#--------------------------------------------------------------------------
+
+
+base_resolution = 1 # Use this to coarsen or refine entire mesh. 
+
+# Basic geometry and bounding polygon
+xleft   = 0
+xright  = 5.448
+ybottom = 0
+ytop    = 3.402
+
+point_sw = [xleft, ybottom]
+point_se = [xright, ybottom]
+point_nw = [xleft, ytop]    
+point_ne = [xright, ytop]
+
+bounding_polygon = [point_se,
+                    point_ne,
+                    point_nw,
+                    point_sw]
+
+
+# Localised refined area for gulleys
+xl = 4.8
+xr = 5.3
+yb = 1.6
+yt = 2.3
+p0 = [xl, yb]
+p1 = [xr, yb]
+p2 = [xr, yt]
+p3 = [xl, yt]
+
+gulleys = [p0, p1, p2, p3]
+
+
+
+# Island area and drawdown region
+island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
+island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
+island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
+island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
+island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
+island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
+island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
+island_7 = [xl-.01, yt]
+
+island = [island_0, island_1, island_2,
+          island_3, island_4, island_5,
+          island_6, island_7]
+
+
+
+# Region spanning half right hand side of domain just inside boundary
+rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
+rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
+rhs_se = [xright-0.1, ybottom+0.2]
+rhs_ne = [xright-0.1, ytop-0.2]        
+
+rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
+
+
+# Interior regions and creation of mesh
+interior_regions = [[rhs_region, 0.0005],
+                    [island, 0.0003*base_resolution],
+                    [gulleys, 0.00003*base_resolution]]    
+
+
 
 
@@ -69,4 +145,6 @@
 
 
+
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
@@ -77,74 +155,4 @@
 
 
-    #--------------------------------------------------------------------------
-    # Create the triangular mesh based on overall clipping polygon with a
-    # tagged
-    # boundary and interior regions defined in project.py along with
-    # resolutions (maximal area of per triangle) for each polygon
-    #--------------------------------------------------------------------------
-
-
-    base_resolution = 1 # Use this to coarsen or refine entire mesh. 
-
-    # Basic geometry and bounding polygon
-    xleft   = 0
-    xright  = 5.448
-    ybottom = 0
-    ytop    = 3.402
-
-    point_sw = [xleft, ybottom]
-    point_se = [xright, ybottom]
-    point_nw = [xleft, ytop]    
-    point_ne = [xright, ytop]
-
-    bounding_polygon = [point_se,
-                        point_ne,
-                        point_nw,
-                        point_sw]
-    
-
-    # Localised refined area for gulleys
-    xl = 4.8
-    xr = 5.3
-    yb = 1.6
-    yt = 2.3
-    p0 = [xl, yb]
-    p1 = [xr, yb]
-    p2 = [xr, yt]
-    p3 = [xl, yt]
-    
-    gulleys = [p0, p1, p2, p3]
-    
-
-
-    # Island area and drawdown region
-    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
-    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
-    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
-    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
-    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
-    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
-    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
-    island_7 = [xl-.01, yt]
- 
-    island = [island_0, island_1, island_2,
-              island_3, island_4, island_5,
-              island_6, island_7]
-
-
-
-    # Region spanning half right hand side of domain just inside boundary
-    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
-    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
-    rhs_se = [xright-0.1, ybottom+0.2]
-    rhs_ne = [xright-0.1, ytop-0.2]        
-
-    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
-
-    
-    # Interior regions and creation of mesh
-    interior_regions = [[rhs_region, 0.0005],
-                        [island, 0.0003*base_resolution],
-                        [gulleys, 0.00003*base_resolution]]    
 
     meshname = project.mesh_filename + '.msh'
