Index: cumentation/user_manual/examples/runup_convergence.py
===================================================================
--- /documentation/user_manual/examples/runup_convergence.py	(revision 3505)
+++ 	(revision )
@@ -1,210 +1,0 @@
-"""Simple water flow example using ANUGA
-
-Water driven up a linear slope and time varying boundary,
-similar to a beach environment.
-
-The study area is discretised as a regular triangular grid 100m x 100m
-"""
-
-
-#------------------------------------------------------------------------------
-# Import necessary modules
-#------------------------------------------------------------------------------
-
-import sys
-
-from pmesh.mesh_interface import create_mesh_from_regions
-from pyvolution.mesh_factory import rectangular_cross
-from pyvolution.shallow_water import Domain
-from pyvolution.shallow_water import Reflective_boundary
-from pyvolution.shallow_water import Dirichlet_boundary
-from pyvolution.shallow_water import Time_boundary
-from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from pyvolution.util import file_function
-from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis
-ion()
-
-
-#------------------------------------------------------------------------------
-# Model constants
-
-slope = -0.02     # 1:50 Slope
-highest_point = 3 # Highest elevation (m)
-sea_level = -2    # Mean sea level
-wave_height = 3   # Solitary wave height H
-simulation_name = 'runup_convergence'   
-
-
-# Basin dimensions (m)
-west = 0          # left boundary
-east = 500        # right boundary 
-south = 0         # lower boundary
-north = 100       # upper bourdary
-
-
-#------------------------------------------------------------------------------
-# Setup computational domain all units in meters
-#------------------------------------------------------------------------------
-
-# Structured mesh
-dx = 10           # Resolution: Lenght of subdivisions on x axis (length)
-dy = 10           # Resolution: Lenght of subdivisions on y axis (width)
-
-length = east-west
-width = north-south
-points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width,
-                                               origin = (west, south)) 
-
-domain = Domain(points, vertices, boundary) # Create domain
-
-
-
-# Unstructured mesh
-#polygon = [[east,north],[west,north],[west,south],[east,south]]
-#meshname = simulation_name + '.msh'
-#create_mesh_from_regions(polygon,
-#                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
-#                         maximum_triangle_area=dx*dy/2,
-#                         filename=meshname)
-#                         #interior_regions=interior_regions)
-#domain = Domain(meshname, use_cache=True, verbose = True)
-
-
-
-domain.set_name(simulation_name)
-
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
-
-def topography(x,y):
-    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
-
-domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0.0 )        # Constant friction 
-domain.set_quantity('stage', sea_level)      # Constant initial stage
-
-
-#------------------------------------------------------------------------------
-# Setup boundary conditions
-#------------------------------------------------------------------------------
-
-from math import sin, pi, cosh, sqrt
-Br = Reflective_boundary(domain)      # Solid reflective wall
-Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
-
-def waveform(t): 
-    return wave_height/cosh(t-5)**2
-
-Bw = Time_boundary(domain=domain,     # Time dependent boundary  
-                   f=lambda t: [waveform(t), 0.0, 0.0])
-
-# Time dependent boundary for stage, where momentum is set automatically
-Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
-
-
-# Associate boundary tags with boundary objects
-domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
-
-
-#------------------------------------------------------------------------------
-# Evolve system through time
-#------------------------------------------------------------------------------
-
-stagestep = []
-for t in domain.evolve(yieldstep = 1, finaltime = 100):
-    domain.write_time()
-
-
-#-----------------------------------------------------------------------------
-# Interrogate solution
-#-----------------------------------------------------------------------------
-
-
-# Define line of gauges through center of domain
-def gauge_line(west,east,north,south):
-    from Numeric import arange
-    x_vector = arange(west, (east-west)/2, 10) # Gauges every 1 meter from west to midpt
-    y = (north+south)/2.
-
-    gauges = []
-    for x in x_vector:
-        gauges.append([x,y])
-        
-    return gauges, x_vector
-
-gauges, x_vector = gauge_line(west,east,north,south)
-
-# Obtain interpolated timeseries at gauges
-f = file_function(domain.get_name()+'.sww',
-                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
-                  interpolation_points = gauges,
-                  verbose = True,
-                  use_cache = True)
-
-
-# Find runup distance from western boundary through a linear search
-max_stage = []
-min_stage = []
-runup_point = west
-coastline = east        
-for k, g in enumerate(gauges):
-    z = f(0, point_id=k)[1] # Elevation
-
-    min_w = sys.maxint
-    max_w = -min_w
-    for i, t in enumerate(f.get_time()):
-        w = f(t, point_id = k)[0]
-        if w > max_w: max_w = w
-        if w < min_w: min_w = w        
-
-    if max_w-z <= 0.01:  # Find first gauge where max depth > eps (runup)
-        runup_point = g[0]
-
-    if min_w-z <= 0.01:  # Find first gauge where min depth > eps (coastline)
-        coastline = g[0]        
-        
-    max_stage.append(max_w)
-    min_stage.append(min_w)    
-
-
-# Print
-
-print 'wave height [m]:                    ', wave_height
-runup_height = topography(runup_point, (north+south)/2.)
-print 'run up height [m]:                  ', runup_height 
-
-runup_distance = runup_point-coastline
-print 'run up distance from coastline [m]: ', runup_distance
-
-print 'Coastline (meters form west):       ', coastline
-
-
-
-# Store
-filename = 'maxrunup'+str(wave_height)+'.csv'
-fid = open(filename,'w')    
-s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
-fid.write(s)
-
-s = '%.2f,%.2f,%.2f\n' %(wave_height, runup_distance, runup_height)
-fid.write(s)
-
-fid.close()
-
-# Plot
-figure(1)
-plot(x_vector, max_stage, 'g+',
-     x_vector, min_stage, 'b+',     
-     x_vector, topography(x_vector,(north+south)/2.), 'r-')
-xlabel('x')
-ylabel('stage')
-title('Maximum stage for gauge line')
-#axis([33000, 47000, -1000, 3000])
-savefig('max_stage')
-
-close('all')
-    
-
-
Index: /inundation/examples/runup_convergence.py
===================================================================
--- /inundation/examples/runup_convergence.py	(revision 3506)
+++ /inundation/examples/runup_convergence.py	(revision 3506)
@@ -0,0 +1,210 @@
+"""Simple water flow example using ANUGA
+
+Water driven up a linear slope and time varying boundary,
+similar to a beach environment.
+
+The study area is discretised as a regular triangular grid 100m x 100m
+"""
+
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+
+from pmesh.mesh_interface import create_mesh_from_regions
+from pyvolution.mesh_factory import rectangular_cross
+from pyvolution.shallow_water import Domain
+from pyvolution.shallow_water import Reflective_boundary
+from pyvolution.shallow_water import Dirichlet_boundary
+from pyvolution.shallow_water import Time_boundary
+from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from pyvolution.util import file_function
+from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis
+ion()
+
+
+#------------------------------------------------------------------------------
+# Model constants
+
+slope = -0.02     # 1:50 Slope
+highest_point = 3 # Highest elevation (m)
+sea_level = -2    # Mean sea level
+wave_height = 3   # Solitary wave height H
+simulation_name = 'runup_convergence'   
+
+
+# Basin dimensions (m)
+west = 0          # left boundary
+east = 500        # right boundary 
+south = 0         # lower boundary
+north = 100       # upper bourdary
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain all units in meters
+#------------------------------------------------------------------------------
+
+# Structured mesh
+dx = 10           # Resolution: Lenght of subdivisions on x axis (length)
+dy = 10           # Resolution: Lenght of subdivisions on y axis (width)
+
+length = east-west
+width = north-south
+points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width,
+                                               origin = (west, south)) 
+
+domain = Domain(points, vertices, boundary) # Create domain
+
+
+
+# Unstructured mesh
+#polygon = [[east,north],[west,north],[west,south],[east,south]]
+#meshname = simulation_name + '.msh'
+#create_mesh_from_regions(polygon,
+#                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
+#                         maximum_triangle_area=dx*dy/2,
+#                         filename=meshname)
+#                         #interior_regions=interior_regions)
+#domain = Domain(meshname, use_cache=True, verbose = True)
+
+
+
+domain.set_name(simulation_name)
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+
+def topography(x,y):
+    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
+
+domain.set_quantity('elevation', topography) # Use function for elevation
+domain.set_quantity('friction', 0.0 )        # Constant friction 
+domain.set_quantity('stage', sea_level)      # Constant initial stage
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+from math import sin, pi, cosh, sqrt
+Br = Reflective_boundary(domain)      # Solid reflective wall
+Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
+
+def waveform(t): 
+    return wave_height/cosh(t-5)**2
+
+Bw = Time_boundary(domain=domain,     # Time dependent boundary  
+                   f=lambda t: [waveform(t), 0.0, 0.0])
+
+# Time dependent boundary for stage, where momentum is set automatically
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+stagestep = []
+for t in domain.evolve(yieldstep = 1, finaltime = 100):
+    domain.write_time()
+
+
+#-----------------------------------------------------------------------------
+# Interrogate solution
+#-----------------------------------------------------------------------------
+
+
+# Define line of gauges through center of domain
+def gauge_line(west,east,north,south):
+    from Numeric import arange
+    x_vector = arange(west, (east-west)/2, 10) # Gauges every 1 meter from west to midpt
+    y = (north+south)/2.
+
+    gauges = []
+    for x in x_vector:
+        gauges.append([x,y])
+        
+    return gauges, x_vector
+
+gauges, x_vector = gauge_line(west,east,north,south)
+
+# Obtain interpolated timeseries at gauges
+f = file_function(domain.get_name()+'.sww',
+                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
+                  interpolation_points = gauges,
+                  verbose = True,
+                  use_cache = True)
+
+
+# Find runup distance from western boundary through a linear search
+max_stage = []
+min_stage = []
+runup_point = west
+coastline = east        
+for k, g in enumerate(gauges):
+    z = f(0, point_id=k)[1] # Elevation
+
+    min_w = sys.maxint
+    max_w = -min_w
+    for i, t in enumerate(f.get_time()):
+        w = f(t, point_id = k)[0]
+        if w > max_w: max_w = w
+        if w < min_w: min_w = w        
+
+    if max_w-z <= 0.01:  # Find first gauge where max depth > eps (runup)
+        runup_point = g[0]
+
+    if min_w-z <= 0.01:  # Find first gauge where min depth > eps (coastline)
+        coastline = g[0]        
+        
+    max_stage.append(max_w)
+    min_stage.append(min_w)    
+
+
+# Print
+
+print 'wave height [m]:                    ', wave_height
+runup_height = topography(runup_point, (north+south)/2.)
+print 'run up height [m]:                  ', runup_height 
+
+runup_distance = runup_point-coastline
+print 'run up distance from coastline [m]: ', runup_distance
+
+print 'Coastline (meters form west):       ', coastline
+
+
+
+# Store
+filename = 'maxrunup'+str(wave_height)+'.csv'
+fid = open(filename,'w')    
+s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
+fid.write(s)
+
+s = '%.2f,%.2f,%.2f\n' %(wave_height, runup_distance, runup_height)
+fid.write(s)
+
+fid.close()
+
+# Plot
+figure(1)
+plot(x_vector, max_stage, 'g+',
+     x_vector, min_stage, 'b+',     
+     x_vector, topography(x_vector,(north+south)/2.), 'r-')
+xlabel('x')
+ylabel('stage')
+title('Maximum stage for gauge line')
+#axis([33000, 47000, -1000, 3000])
+savefig('max_stage')
+
+close('all')
+    
+
+
