Index: /anuga_validation/okushiri_2005/create_okushiri.py
===================================================================
--- /anuga_validation/okushiri_2005/create_okushiri.py	(revision 3853)
+++ /anuga_validation/okushiri_2005/create_okushiri.py	(revision 3854)
@@ -1,8 +1,10 @@
 """Create mesh and time boundary for the Okushiri island validation
 """
+
 
 from Numeric import array, zeros, Float, allclose
 
 from anuga.pmesh.mesh import *
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 
@@ -84,8 +86,14 @@
 
 
-    # Create Mesh
+    #--------------------------------------------------------------------------
+    # 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
+    #--------------------------------------------------------------------------
 
-    #Basic geometry
-    
+    base_resolution = 10
+
+    # Basic geometry
     xleft   = 0
     xright  = 5.448
@@ -93,5 +101,5 @@
     ytop    = 3.402
 
-    #Outline 
+    # Outline 
     point_sw = [xleft, ybottom]
     point_se = [xright, ybottom]
@@ -99,50 +107,9 @@
     point_ne = [xright, ytop]
 
-    #Midway points (left)
+    # Midway points (left)
     point_mtop = [xleft + (xright-xleft)/3+1, ytop]
     point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
 
-
-    geo = Geo_reference(xllcorner = xleft,
-                        yllcorner = ybottom)
-    
-                        
-    print "***********************"
-    print "geo ref", geo
-    print "***********************"
-    
-    m = Mesh(geo_reference=geo)
-
-    #Boundary
-    dict = {}
-    dict['points'] = [point_se,   #se
-                      point_ne,
-                      point_mtop,
-                      point_nw,
-                      point_sw,
-                      point_mbottom]
-
-    
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
-                        [4,5], [5,0],  #The outer border
-                        [2,5]]         #Separator
-    
-    dict['segment_tags'] = ['wall',
-                            'wall',
-                            'wall',
-                            'wave',
-                            'wall',
-                            'wall',
-                            '']           #Interior
-
-        
-    m.addVertsSegs(dict)
-
-
-
-
-
     #Localised refined area for gulleys
-    dict = {}
     xl = 4.8
     xr = 5.3
@@ -153,9 +120,5 @@
     p2 = [xr, yt]
     p3 = [xl, yt]
-    
-    dict['points'] = [p0, p1, p2, p3]
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
-    dict['segment_tags'] = ['', '', '', '']
-    m.addVertsSegs(dict)    
+    gulleys = [p0, p1, p2, p3]
 
     #Island area
@@ -164,64 +127,41 @@
     island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
     island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
-    island_4 = [xleft + 2*(xright-xleft)/3+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.2]
     island_6 = [xl-.01, yb]  #OK
     island_7 = [xl-.01, yt]  #OK      
+
+    island = [island_0, island_1, island_2,
+              island_3, island_4, island_5,
+              island_6, island_7]
+
+
     
 
-    dict['points'] = [island_0, island_1, island_2,
-                      island_3, island_4, island_5,
-                      #p0, p3]                      
-                      island_6, island_7]
+    bounding_polygon = [point_se,
+                        point_ne,
+                        point_mtop,
+                        point_nw,
+                        point_sw,
+                        point_mbottom]
+    
 
-                      
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
-                        [4,5], [5,6], [6,7], [7,0]]
-                        
-    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
+    interior_regions = [[gulleys, 0.00002*base_resolution],
+                        [island, 0.00007*base_resolution]]
+    
 
 
-    m.addVertsSegs(dict)    
+    
 
-    
-#
-    
-    base_resolution = 1
+    print 'start create mesh from regions'
 
-    ocean = m.addRegionEN(xleft+.1, ybottom+.1)
-    ocean.setMaxArea(0.1*base_resolution)
+    meshname = project.mesh_filename + '.msh'
+    m = create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
+                                                'wave': [3]},
+                                 maximum_triangle_area=0.1*base_resolution,
+                                 interior_regions=interior_regions,
+                                 filename=project.mesh_filename,
+                                 use_cache=True,
+                                 verbose=True)
 
-    mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
-    mid.setMaxArea(0.0005*base_resolution)
-  
-
-    inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
-    inner.setMaxArea(0.00007*base_resolution)
-
-
-    gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
-    gulleys.setMaxArea(0.00002*base_resolution)
-
-
-    # From r 1709 11 August 2005
-    #base_resolution = 100
-    #
-    #ocean = m.addRegionEN(xleft+.1, ybottom+.1)
-    #ocean.setMaxArea(0.01*base_resolution)
-    #
-    #mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
-    #mid.setMaxArea(0.001*base_resolution)
-
-    #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
-    #inner.setMaxArea(0.0001*base_resolution)
-
-
-    #gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
-    #gulleys.setMaxArea(0.00001*base_resolution)        
-    
-    
-    m.generateMesh('pzq28.0za1000000a')
-
-    import project
-    m.export_mesh_file(project.mesh_filename)
-    
Index: /anuga_validation/okushiri_2005/create_okushiri_original.py
===================================================================
--- /anuga_validation/okushiri_2005/create_okushiri_original.py	(revision 3854)
+++ /anuga_validation/okushiri_2005/create_okushiri_original.py	(revision 3854)
@@ -0,0 +1,227 @@
+"""Create mesh and time boundary for the Okushiri island validation
+"""
+
+from Numeric import array, zeros, Float, allclose
+
+from anuga.pmesh.mesh import *
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+
+import project
+
+
+def prepare_timeboundary(filename):
+    """Converting benchmark 2 time series to NetCDF tms file.
+    This is a 'throw-away' code taylor made for files like
+    'Benchmark_2_input.txt' from the LWRU2 benchmark
+    """
+
+    textversion = filename[:-4] + '.txt'
+
+    print 'Preparing time boundary from %s' %textversion
+    from Numeric import array
+
+    fid = open(textversion)
+
+    #Skip first line
+    line = fid.readline()
+
+    #Read remaining lines
+    lines = fid.readlines()
+    fid.close()
+
+
+    N = len(lines)
+    T = zeros(N, Float)  #Time
+    Q = zeros(N, Float)  #Values
+
+    for i, line in enumerate(lines):
+        fields = line.split()
+
+        T[i] = float(fields[0])
+        Q[i] = float(fields[1])
+
+
+    #Create tms file
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    print 'Writing to', filename
+    fid = NetCDFFile(filename, 'w')
+
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'Input wave for Benchmark 2'
+    fid.starttime = 0.0
+    fid.createDimension('number_of_timesteps', len(T))
+    fid.createVariable('time', Float, ('number_of_timesteps',))
+    fid.variables['time'][:] = T
+
+    fid.createVariable('stage', Float, ('number_of_timesteps',))
+    fid.variables['stage'][:] = Q[:]
+
+    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
+    fid.variables['xmomentum'][:] = 0.0
+
+    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
+    fid.variables['ymomentum'][:] = 0.0
+
+    fid.close()
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+
+    #Preparing points
+    #(Only when using original Benchmark_2_Bathymetry.txt file)
+    # FIXME: Replace by using geospatial data
+    #
+    #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
+    #xya2pts(project.bathymetry_filename, verbose = True,
+    #        z_func = lambda z: -z)
+
+
+    # Prepare time boundary
+    prepare_timeboundary(project.boundary_filename)
+
+
+    # Create Mesh
+
+    #Basic geometry
+    
+    xleft   = 0
+    xright  = 5.448
+    ybottom = 0
+    ytop    = 3.402
+
+    #Outline 
+    point_sw = [xleft, ybottom]
+    point_se = [xright, ybottom]
+    point_nw = [xleft, ytop]    
+    point_ne = [xright, ytop]
+
+    #Midway points (left)
+    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
+    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
+
+
+    geo = Geo_reference(xllcorner = xleft,
+                        yllcorner = ybottom)
+    
+                        
+    print "***********************"
+    print "geo ref", geo
+    print "***********************"
+    
+    m = Mesh(geo_reference=geo)
+
+    #Boundary
+    dict = {}
+    dict['points'] = [point_se,   #se
+                      point_ne,
+                      point_mtop,
+                      point_nw,
+                      point_sw,
+                      point_mbottom]
+
+    
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
+                        [4,5], [5,0],  #The outer border
+                        [2,5]]         #Separator
+    
+    dict['segment_tags'] = ['wall',
+                            'wall',
+                            'wall',
+                            'wave',
+                            'wall',
+                            'wall',
+                            '']           #Interior
+
+        
+    m.addVertsSegs(dict)
+
+
+
+
+
+    #Localised refined area for gulleys
+    dict = {}
+    xl = 4.8
+    xr = 5.3
+    yb = 1.6
+    yt = 2.3
+    p0 = [xl, yb]
+    p1 = [xr, yb]
+    p2 = [xr, yt]
+    p3 = [xl, yt]
+    
+    dict['points'] = [p0, p1, p2, p3]
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
+    dict['segment_tags'] = ['', '', '', '']
+    m.addVertsSegs(dict)    
+
+    #Island area
+    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.3, ybottom + 2*(ytop-ybottom)/3-0.3]
+    island_3 = [xleft + (xright-xleft)/2+0.3, 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.2]
+    island_6 = [xl-.01, yb]  #OK
+    island_7 = [xl-.01, yt]  #OK      
+    
+
+    dict['points'] = [island_0, island_1, island_2,
+                      island_3, island_4, island_5,
+                      #p0, p3]                      
+                      island_6, island_7]
+
+                      
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
+                        [4,5], [5,6], [6,7], [7,0]]
+                        
+    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
+
+
+    m.addVertsSegs(dict)    
+
+    
+#
+    
+    base_resolution = 1
+
+    ocean = m.addRegionEN(xleft+.1, ybottom+.1)
+    ocean.setMaxArea(0.1*base_resolution)
+
+    mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
+    mid.setMaxArea(0.0005*base_resolution)
+  
+
+    inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
+    inner.setMaxArea(0.00007*base_resolution)
+
+
+    gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
+    gulleys.setMaxArea(0.00002*base_resolution)
+
+
+    # From r 1709 11 August 2005
+    #base_resolution = 100
+    #
+    #ocean = m.addRegionEN(xleft+.1, ybottom+.1)
+    #ocean.setMaxArea(0.01*base_resolution)
+    #
+    #mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
+    #mid.setMaxArea(0.001*base_resolution)
+
+    #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
+    #inner.setMaxArea(0.0001*base_resolution)
+
+
+    #gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
+    #gulleys.setMaxArea(0.00001*base_resolution)        
+    
+    
+    m.generateMesh('pzq28.0za1000000a')
+
+    import project
+    m.export_mesh_file(project.mesh_filename)
+    
Index: /anuga_validation/okushiri_2005/extract_timeseries.py
===================================================================
--- /anuga_validation/okushiri_2005/extract_timeseries.py	(revision 3853)
+++ /anuga_validation/okushiri_2005/extract_timeseries.py	(revision 3854)
@@ -21,4 +21,5 @@
     plotting = True
 
+plotting = False
 
 #-------------------------
@@ -42,4 +43,10 @@
                         'ch7': 1.121816242516758e-04,
                         'ch9': 1.249543278366778e-04}
+
+# old limiters
+expected_covariances = {'Boundary': 5.288601162783020386e-05,
+                        'ch5': 1.167001054284431472e-04,
+                        'ch7': 1.121474766904651861e-04,
+                        'ch9': 1.249244820847215335e-04}
 
 
@@ -81,5 +88,5 @@
 # Read and interpolate model output
 #--------------------------------------------------
-#f = file_function('okushiri_new_limiters.sww',
+#f = file_function('okushiri_new_limiters.sww',   #The best so far
 #f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
 f = file_function(project.output_filename,
@@ -111,8 +118,9 @@
         print 'Result is better than expected by: %.18e'\
               %(res-expected_covariances[name])
-        print 'Result = %.18e' %res
         print 'Expect = %.18e' %expected_covariances[name]    
     elif res > expected_covariances[name]+eps:
-        print 'FAIL: %.18e' %res
+        print 'FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_covariances[name])
+        print 'Expect = %.18e' %expected_covariances[name]
     else:
         pass
Index: /anuga_validation/okushiri_2005/project.py
===================================================================
--- /anuga_validation/okushiri_2005/project.py	(revision 3853)
+++ /anuga_validation/okushiri_2005/project.py	(revision 3854)
@@ -15,5 +15,6 @@
 
 # Model output
-output_filename = 'okushiri_old_limiters.sww'
+#output_filename = 'okushiri_old_limiters.sww'
+output_filename = 'okushiri_nl_new_mesh.sww'
 
 
Index: /anuga_validation/okushiri_2005/run_okushiri.py
===================================================================
--- /anuga_validation/okushiri_2005/run_okushiri.py	(revision 3853)
+++ /anuga_validation/okushiri_2005/run_okushiri.py	(revision 3854)
@@ -53,5 +53,5 @@
 domain.set_name(project.output_filename)  # Name of output sww file
 domain.set_default_order(2)               # Apply second order scheme 
-domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
+#domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
 domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
 domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
