Index: /inundation/examples/beach.py
===================================================================
--- /inundation/examples/beach.py	(revision 2619)
+++ /inundation/examples/beach.py	(revision 2620)
@@ -29,5 +29,5 @@
 print 'Creating domain from %s.tsh' %name
 #domain = pmesh_to_domain_instance(name + '.tsh', Domain)
-domain = Domain(mesh_file= name + '.tsh')
+domain = Domain(mesh_filename = name + '.tsh')
 
 domain.store = True
Index: /inundation/examples/island.py
===================================================================
--- /inundation/examples/island.py	(revision 2619)
+++ /inundation/examples/island.py	(revision 2620)
@@ -17,4 +17,5 @@
 from pyvolution.mesh_factory import rectangular
 from pyvolution.shallow_water import Domain, Reflective_boundary
+from pmesh.mesh_interface import create_mesh_from_regions
 from utilities.polygon import Polygon_function, read_polygon
 
@@ -28,6 +29,23 @@
 points, vertices, boundary = rectangular(N, N, 100, 100)
 
+
+create_mesh_from_regions( [[0,0], [100,0], [100,100], [0,100]],
+                          boundary_tags = {'bottom': [0],
+                                           'right': [1],
+                                           'top': [2],
+                                           'left': [3]},
+                          maximum_triangle_area = 25,
+                          filename = 'island.msh',
+                          interior_regions=[ ([[55,25], [75,25],
+                                               [75,75], [55,75]], 3)]
+                          )
+                          
+                                                    
+                          
+                           
+
 #Create shallow water domain
-domain = Domain(points, vertices, boundary)
+#domain = Domain(points, vertices, boundary)
+domain = Domain(mesh_filename = 'island.msh')
 domain.smooth = False
 domain.set_name('island')
@@ -55,5 +73,5 @@
     return z
 
-domain.set_quantity('friction', 0.1)
+domain.set_quantity('friction', 1.0)
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
@@ -75,5 +93,5 @@
 
 import time
-for t in domain.evolve(yieldstep = 1, finaltime = 100):
+for t in domain.evolve(yieldstep = 0.5, finaltime = 100):
     domain.write_time()
 
Index: /inundation/examples/island_timeseries.py
===================================================================
--- /inundation/examples/island_timeseries.py	(revision 2620)
+++ /inundation/examples/island_timeseries.py	(revision 2620)
@@ -0,0 +1,195 @@
+"""Read in sww file, interpolate at specified locations and plot time series
+
+"""
+
+from pyvolution.util import file_function
+from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
+from pylab import *
+from matplotlib.ticker import MultipleLocator, FormatStrFormatter
+
+
+swwfile = 'island.sww' 
+#gauges = [[0,0.5], [0.2,0.5], [0.4,0.5], [0.6,0.5], [0.8,0.5], [1,0.5]]
+#gauges = [[20, 50], [30, 50], [40, 50], [50, 50]]
+gauges = [[42, 48], [58, 48]]     
+
+
+#Read model output
+quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+f = file_function(swwfile,
+                  quantities = quantities,
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+               
+from math import sqrt, atan, degrees
+from Numeric import ones
+N = len(gauges)
+for k, g in enumerate(gauges):
+    if k%((N+10)/10)==0: # diagnostics - print 10 lines
+        print 'Doing row %d of %d' %(k, N)
+
+    model_time = []
+    stages = []
+    elevations = []
+    momenta = []
+    velocity = []
+    xmom = []
+    ymom = []
+    bearings = []
+    depths = []
+
+    max_depth = 0
+    max_momentum = 0
+    max_velocity = 0
+
+    due_east = 90.0*ones([len(f.T)])
+    due_west = 270.0*ones([len(f.T)])
+    maxT = max(f.T)
+    tstep = maxT/(len(f.T)-1)
+
+
+    myloc = 'G%d' %k #locations[k]
+
+    for i, t in enumerate(f.T): # T is a list of times
+        #if tmin < t < tmax:
+        w = f(t, point_id = k)[0]
+        z = f(t, point_id = k)[1]
+        uh = f(t, point_id = k)[2]
+        vh = f(t, point_id = k)[3]
+        depth = w-z
+
+        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
+        vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
+        angle = degrees(atan(vh/(uh+1.e-15)))
+        if (0 < angle < 90.0):
+            if vh > 0:
+                bearing = 90.0 - abs(angle)
+            if vh < 0:
+                bearing = 270.0 - abs(angle)
+        if (-90 < angle < 0):
+            if vh < 0:
+                bearing = 90.0 - abs(angle)
+            if vh > 0:
+                bearing = 270.0 - abs(angle)
+        if angle == 0:
+            bearing = 0.0
+                
+        model_time.append(t)        
+        stages.append(w)
+        elevations.append(z)  #Should be constant over time
+        momenta.append(m)
+        velocity.append(vel)
+        xmom.append(uh)
+        ymom.append(vh)
+        bearings.append(bearing)
+        depths.append(depth)
+
+        if w-z > max_depth:
+            max_depth = w-z
+        if m > max_momentum:            
+            max_momentum = m
+        if vel > max_velocity:
+            max_velocity = vel
+
+    #Plot only those gauges that have been inundated by more than a threshold
+    #if max_depth < 0.2:
+    #    print 'Skipping gauge %d' %k
+    #    continue
+
+    ion()
+    hold(False)
+
+    if elevations[0] < -10:
+        plot(model_time, stages, '-b')
+    else:    
+        plot(model_time, stages, '-b',
+             model_time, elevations, '-k')
+    axis([0, 100, z-0.01, z+0.005])
+        
+    #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
+    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
+    title(name)
+
+    title('%s (stage)' %name)
+    xlabel('time [s]')
+    ylabel('elevation [m]')    
+    legend(('Stage', 'Bed = %.1f' %elevations[0]),
+           shadow=True,
+           loc='upper right')
+    #savefig('Gauge_%d_stage' %k)
+    savefig('Gauge_%s_stage' %myloc)
+
+    raw_input('Next')
+
+    """
+    #Momentum plot
+    ion()
+    hold(False)
+    plot(model_time, momenta, '-r')
+    title(name)
+
+    title('%s (momentum)' %name)
+    xlabel('time [s]')
+    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')    
+    #savefig('Gauge_%d_momentum' %k)
+    savefig('Gauge_%s_momentum' %myloc)
+    
+    raw_input('Next')
+    """
+
+    """
+    #Bearing plot
+    ion()
+    hold(False)
+    ax = plot(model_time, bearings, '-b', model_time, due_west, '-.b',
+         model_time, due_east, '-.b')
+    title(name)
+    ax = axis([0, maxT, 0, 360])
+    text(maxT+tstep, 90, 'East')
+    text(maxT+tstep, 270, 'West')
+    #majorLocator = MultipleLocator(3600)
+    #print 'major', majorLocator[1]
+    #ax.xaxis.set_major_locator(majorLocator) #'yticklabels', range(30,390,30))
+    # set(labels,color='g',rotation=45)
+
+    title('%s (bearing)' %name)
+    xlabel('time [s]')
+    ylabel(' atan(vh/uh) [degrees from North]')    
+    #savefig('Gauge_%d_bearing' %k)
+    savefig('Gauge_%s_bearing' %myloc)
+    
+    raw_input('Next')
+    """
+
+    """
+    #Speed plot
+    ion()
+    hold(False)
+    plot(model_time, velocity, '-r')
+    title(name)
+
+    title('%s (velocity)' %name)
+    xlabel('time [s]')
+    ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')    
+    #savefig('Gauge_%d_speed' %k)
+    savefig('Gauge_%s_speed' %myloc)
+    
+    raw_input('Next')
+    """
+    
+
+    whichone = '_%s' %myloc
+    thisfile = 'island_gauge'+whichone+'.csv'
+    fid = open(thisfile, 'w')
+    s = 'Time, Depth, Momentum, Velocity \n'
+    fid.write(s)
+    for i_t, i_d, i_m, i_vel in zip(model_time, depths, momenta, velocity):
+        s = '%.2f, %.2f, %.2f, %.2f\n' %(i_t, i_d, i_m, i_vel)
+        fid.write(s)
+        
+show()
+
+
Index: /inundation/pyvolution/domain.py
===================================================================
--- /inundation/pyvolution/domain.py	(revision 2619)
+++ /inundation/pyvolution/domain.py	(revision 2620)
@@ -20,4 +20,8 @@
 class Domain(Mesh):
 
+
+    #FIXME (Ole): I would like this interface to automatically work out
+    #whether input is mesh data or a filename (as in set_quantity).
+    #That way the user need not remember the keyword arguments.
     def __init__(self, coordinates=None, vertices=None,
                  boundary=None,
@@ -25,10 +29,10 @@
                  tagged_elements=None, geo_reference=None,
                  use_inscribed_circle=False,
-                 mesh_file=None):
-
-        if mesh_file is not None:
+                 mesh_filename=None):
+
+        if mesh_filename is not None:
             coordinates, vertices, boundary, vertex_quantity_dict \
                                 ,tagged_elements, geo_reference = \
-                                pmesh_to_domain(file_name=mesh_file)
+                                pmesh_to_domain(file_name=mesh_filename)
             # FIXME(DSG-DSG) have to do something with vertex_quantity_dict
 
@@ -109,5 +113,5 @@
         self.already_computed_flux = zeros((N, 3), Int)
 
-        if mesh_file is not None:
+        if mesh_filename is not None:
             # If the mesh file passed any quantity values
             # , initialise with these values. 
Index: /inundation/pyvolution/shallow_water.py
===================================================================
--- /inundation/pyvolution/shallow_water.py	(revision 2619)
+++ /inundation/pyvolution/shallow_water.py	(revision 2620)
@@ -72,5 +72,5 @@
                  tagged_elements = None, geo_reference = None,
                  use_inscribed_circle=False,
-                 mesh_file=None):
+                 mesh_filename=None):
 
         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
@@ -78,5 +78,6 @@
         Generic_Domain.__init__(self, coordinates, vertices, boundary,
                                 conserved_quantities, other_quantities,
-                                tagged_elements, geo_reference, use_inscribed_circle, mesh_file)
+                                tagged_elements, geo_reference,
+                                use_inscribed_circle, mesh_filename)
 
         from config import minimum_allowed_height, maximum_allowed_speed, g
Index: /inundation/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/pyvolution/test_shallow_water.py	(revision 2619)
+++ /inundation/pyvolution/test_shallow_water.py	(revision 2620)
@@ -3329,5 +3329,5 @@
          #domain = pmesh_to_domain_instance(fileName, Domain)
          
-         domain = Domain(mesh_file=fileName)
+         domain = Domain(mesh_filename=fileName)
          
          os.remove(fileName)
